{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Least Angle Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 主要实现LARS算法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Forward-stagewise regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Alg**: <br/>\n",
    "$initialization\\ \\hat{\\mu} = 0$ <br/>\n",
    "$repeat\\ for\\ t = 1, 2, 3, ...:$<br/>\n",
    "$\\quad Cal\\ current\\ correlations\\ vector\\ \\hat{c}$ <br/>\n",
    "$\\quad Pick\\ up\\ greatest\\ current\\ correlations\\ variable$ <br/>\n",
    "$\\quad Take\\ a\\ small\\ step\\ 0 < \\epsilon < |c_j|\\ in\\ the\\ direction\\ of\\ selected\\ variable$ <br/>\n",
    "$\\quad \\hat{\\mu} \\leftarrow \\hat{\\mu} + \\epsilon * sign(c_j) * x_j$ <br/>\n",
    "$until\\ none\\ of\\ the\\ variables\\ have\\ correlation\\ with\\ the\\ residual$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import re\n",
    "import numpy as np\n",
    "from sklearn import linear_model\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def forward_stagewise_reg(X, y, t=5000, epsilon=3.3):\n",
    "    mu_hat = np.zeros_like(y)\n",
    "    w = np.zeros_like(X[0])\n",
    "    \n",
    "    while True:\n",
    "        # cal 𝑐𝑢𝑟𝑟𝑒𝑛𝑡 𝑐𝑜𝑟𝑟𝑒𝑙𝑎𝑡𝑖𝑜𝑛 vector\n",
    "        C = (X.T).dot(y - mu_hat)\n",
    "        \n",
    "        # pick up highest 𝑐𝑢𝑟𝑟𝑒𝑛𝑡 𝑐𝑜𝑟𝑟𝑒𝑙𝑎𝑡𝑖𝑜𝑛 variable\n",
    "        j = np.argmax(np.abs(C))\n",
    "        \n",
    "        # cal step_length\n",
    "        step_length = epsilon if epsilon < abs(C[j]) else abs(C[j])\n",
    "        \n",
    "        # update mu_hat and w\n",
    "        mu_hat += step_length * np.sign(C[j]) * X[:, j]\n",
    "        w[j] +=  step_length * np.sign(C[j])\n",
    "        \n",
    "        if np.linalg.norm(w, ord=1) > t or np.linalg.norm(C, ord=1) < 1e-6:\n",
    "            break\n",
    "    return w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_data_from_file(fpath):\n",
    "    label = []\n",
    "    data = []\n",
    "    with open(fpath, 'r') as f:\n",
    "        label = f.readline().strip().split('\\t')\n",
    "        for line in f:\n",
    "            data.append(re.split(' |\\t', line.strip()))\n",
    "    data = np.array(data, dtype=np.float64)\n",
    "    return data, label"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(442, 11)\n",
      "[[ 59.       2.      32.1    101.     157.      93.2     38.       4.\n",
      "    4.8598  87.     151.    ]\n",
      " [ 48.       1.      21.6     87.     183.     103.2     70.       3.\n",
      "    3.8918  69.      75.    ]\n",
      " [ 72.       2.      30.5     93.     156.      93.6     41.       4.\n",
      "    4.6728  85.     141.    ]\n",
      " [ 24.       1.      25.3     84.     198.     131.4     40.       5.\n",
      "    4.8903  89.     206.    ]\n",
      " [ 50.       1.      23.     101.     192.     125.4     52.       4.\n",
      "    4.2905  80.     135.    ]]\n"
     ]
    }
   ],
   "source": [
    "data, label = read_data_from_file('diabetes.data')\n",
    "print(data.shape)\n",
    "print(data[:5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.03807591  0.05068012  0.06169621  0.02187239 -0.0442235  -0.03482076\n",
      "  -0.04340085 -0.00259226  0.01990749 -0.01764613]\n",
      " [-0.00188202 -0.04464164 -0.05147406 -0.02632753 -0.00844872 -0.01916334\n",
      "   0.07441156 -0.03949338 -0.06833155 -0.09220405]\n",
      " [ 0.08529891  0.05068012  0.04445121 -0.00567042 -0.04559945 -0.03419447\n",
      "  -0.03235593 -0.00259226  0.00286131 -0.02593034]\n",
      " [-0.08906294 -0.04464164 -0.01159501 -0.03665608  0.01219057  0.02499059\n",
      "  -0.03603757  0.03430886  0.02268774 -0.00936191]\n",
      " [ 0.00538306 -0.04464164 -0.03638469  0.02187239  0.00393485  0.01559614\n",
      "   0.00814208 -0.00259226 -0.03198764 -0.04664087]]\n",
      "[ -1.13348416 -77.13348416 -11.13348416  53.86651584 -17.13348416]\n"
     ]
    }
   ],
   "source": [
    "X, y = data[:, :-1], data[:,-1]\n",
    "X = (X - X.mean(axis=0))\n",
    "X = X / np.linalg.norm(X, axis=0)\n",
    "y -= y.mean()\n",
    "print(X[:5])\n",
    "print(y[:5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "OLS:  [ -10.0098663  -239.81564367  519.84592005  324.3846455  -792.17563855\n",
      "  476.73902101  101.04326794  177.06323767  751.27369956   67.62669218]\n"
     ]
    }
   ],
   "source": [
    "print('OLS: ', (np.linalg.inv((X.T).dot(X))).dot(X.T).dot(y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ -10.00986617, -239.81564348,  519.84592049,  324.38464525,\n",
       "       -792.17561537,  476.73900289,  101.04325735,  177.06323414,\n",
       "        751.27369086,   67.62669233])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "forward_stagewise_reg(X ,y, t=3460)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts = np.linspace(0, 3460, 100)\n",
    "coefs = []\n",
    "\n",
    "for t in ts:\n",
    "    w = forward_stagewise_reg(X ,y, t=t)\n",
    "    coefs.append(w)\n",
    "coefs = np.array(coefs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-173.0, 3633.0, -869.3480806781228, 828.4461561688399)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAt8AAAH5CAYAAABZF1JyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xl8VPW9//HXd2aSySSTbRIgJGQBQ1gjCJElgESlXHCtVI3CvV7ootfeor2t91Zba9Ur92dre20tvRVqH5brFaVitZW6VAUlkLhACLIvsmUFkskyk8lktu/vjxlSQFZJcrJ8no/HPGZyznfOeZ8oyWe++Z7vV2mtEUIIIYQQQnQ9k9EBhBBCCCGE6C+k+BZCCCGEEKKbSPEthBBCCCFEN5HiWwghhBBCiG4ixbcQQgghhBDdRIpvIYQQQgghuokU30II0Y8ppbRSKtfgDAuVUhsudp8QQvRGUnwLIcQZKKUOKaVmGZ2jMymlPlBKfdOgc09XSpUqpZqVUk6l1Eal1JVGZBFCCCNZjA4ghBCib1NKJQBrgHuBPwLRwAyg3cBMFq11wKjzCyH6L+n5FkKIi6CUSlZKrVFKHVdKNUZeDzlp/0Kl1AGllEspdVAptSCyPVcp9WGk57deKbXqpPcUKqU+jez7VClVeI7zH1JKPaSU2hk5//NKqZjzZVNKLSFc8C5VSrmVUktPOuwspdS+yHt+o5RS58t8kfIAtNYvaa2DWus2rfXftNafneUan1JKbVBKJZ5h30il1LuR3vM9SqnbT9p3vVJqi1KqRSlVqZR69KR9OZEhNt9QSh0B1p607Z+VUkci1/ijL3mNQghxQaT4FkKIi2MCngeygSygDVgKoJSKA54B5mqt44FCoCLyvv8E/gYkA0OAX0fe4wD+GnlfCvDfwF+VUinnyLAA+AfgMsKF7cPny6a1/hFQAnxHa23XWn/npOPdAFwJjANujxz7rJm/hL1AUCm1Qik1VymVfKZGSimTUup3wOXAbK1182n744B3gZXAQOBO4H+UUmMiTVqBu4Ak4HrgXqXUV087zUxg1EnXCDAdGAFcCzyilBr1Ja9TCCHOS4pvIYS4CFrrBq31q1prj9baBSwhXNCdEALGKqVsWutarfWOyHY/4aI4XWvt1VqfuInwemCf1voFrXVAa/0SsBu48RwxlmqtK7XWzsj577zAbGfzpNa6SWt9BFgHjD9P5ouitW4hXOBq4HfAcaXUX5RSg05qFgW8BDiAG7XWnjMc6gbgkNb6+cj3qhx4Fbg1cp4PtNbbtNahSK/6S3zx+h/VWrdqrdtO2vZYpDd+K7CV8IcQIYToElJ8CyHERVBKxSqllimlDiulWoD1QJJSyqy1bgWKgX8BapVSf1VKjYy89T8ABXyilNqhlPp6ZHs6cPi00xwGMs4Ro/K0tunny3aey6o76bUHsJ8n8ymUUs9GhrK4lVI/PFMbrfUurfVCrfUQYGwk8y9PapIL3Ey4EPadJWc2MFkp1XTiQfivAGmRHJOVUusiw26aCf93SD3tGJV80dmuXwghOp0U30IIcXG+T3iIwmStdQJwVWS7AtBav6O1/gowmHAP9u8i2+u01t/SWqcD9xAeLpEL1BAuKk+WBVSfI0PmaW1rLiQb4Z7nC3aOzKe3+5fIUBa71vq/LuC4u4E/EC7CT9gFLALeUkqNOMtbK4EPtdZJJz3sWut7I/tXAn8BMrXWicCz/P3aO05/vnxCCNGVpPgWQoizi1JKxZz0sADxhMdSN0XGa//kRGOl1CCl1E2RscntgBsIRvbddtKNmY2Ei8Ag8CaQp5Sar5SyKKWKgdGEZwc5m39VSg2JnP+HwIkbIc+aLeIoMOxCL/4cmS9K5CbJ759082cm4aEyH53cLjLk5ofAe0qpy85wqDWEv1f/pJSKijyuPGmMdjzg1Fp7lVKTgPkXm1UIIbqaFN9CCHF2bxIuZk88HiU8VMIG1BMuHt8+qb2JcO9zDeAkPN7425F9VwIfK6XchHtn79daH9RaNxAey/x9oIHwUI8btNb158i1kvCNkAcijyci28+VDeBXwK2RWU2euYDrP2PmC3jf6VzA5MixWiPZthO+5lNorVcAjxOZjeS0fS5gNnAH4e9xHfBTwBpp8m3gcaWUC3iE8LSGQgjRoyit5S9wQgjRWyilDgHf1Fq/Z3QWIYQQF096voUQQgghhOgmUnwLIYQQQgjRTWTYiRBCCCGEEN1Eer6FEEIIIYToJlJ8CyGEEEII0U0sRgfoSqmpqTonJ8foGEIIIYQQoo/bvHlzvdZ6wPna9eniOycnh02bNhkdQwghhBBC9HFKqcMX0k6GnQghhBBCCNFNDC2+lVL/ppTaoZTarpR6KbJ881Cl1MdKqX1KqVVKqehIW2vk6/2R/TlGZhdCCCGEEOJiGVZ8K6UygPuAAq31WMBMeMngnwJPa62HA43ANyJv+QbQqLXOBZ6OtBNCCCGEEKLXMHrMtwWwKaX8QCxQC1wDzI/sXwE8CvwWuDnyGmA1sFQppfRFTlTu9/upqqrC6/VeevoeKCYmhiFDhhAVFWV0FCGEEEIIcRrDim+tdbVS6ufAEaAN+BuwGWjSWgcizaqAjMjrDKAy8t6AUqoZSAHqL+a8VVVVxMfHk5OTg1KqE66k59Ba09DQQFVVFUOHDjU6jhBCCCGEOI2Rw06SCfdmDwXSgThg7hmanujZPlOl/IVeb6XU3UqpTUqpTcePH//CG7xeLykpKX2u8AZQSpGSktJne/WFEEIIIXo7I2+4nAUc1Fof11r7gT8BhUCSUupEj/wQoCbyugrIBIjsTwScpx9Ua71ca12gtS4YMODMUy32xcL7hL58bUIIIYQQvZ2RxfcRYIpSKlaFK8ZrgZ3AOuDWSJt/Bv4cef2XyNdE9q+92PHePc1rr72GUordu3d3bNu3bx833HADl112GRMnTuTqq69m/fr1APzhD39gwIABjB8/vuOxc+dOo+ILIYQQQoiLZFjxrbX+mPCNk+XAtkiW5cAPgO8ppfYTHtP9+8hbfg+kRLZ/D3iw20N3spdeeonp06fz8ssvA+EhMddffz133303n3/+OZs3b+bXv/41Bw4c6HhPcXExFRUVHY/Ro0cbFV8IIYQQQlwkQ2c70Vr/BPjJaZsPAJPO0NYL3NYdubqD2+1m48aNrFu3jptuuolHH32UF198kalTp3LTTTd1tBs7dixjx441MKkQQgghhOgsRk81aKjH3tjBzpqWTj3m6PQEfnLjmPO2e/3115kzZw55eXk4HA7Ky8vZsWMHEyZMOOf7Vq1axYYNGzq+Lisrw2azXXJuIYQQQgjR9WR5eYO89NJL3HHHHQDccccdvPTSS19oc8sttzB27FjmzZvXse30YSdSeAshhBBC9B79uuf7Qnqou0JDQwNr165l+/btKKUIBoMopfjJT37ScXMlhG/I3LRpEw888IAhOYUQQgghROeSnm8DrF69mrvuuovDhw9z6NAhKisrGTp0KHl5eWzcuJG//OUvHW09Ho+BSYUQQgghRGeS4tsAL730Erfccssp2772ta+xcuVK1qxZw7PPPsuwYcOYOnUqTzzxBA8//HBHu1WrVp0y1WBpaWl3xxdCCCGEEF+S6uVTZZ9TQUGB3rRp0ynbdu3axahRowxK1D36wzUKIYQQQvQkSqnNWuuC87WTnm8hhBBCCNHraa3xttcRCLiMjnJO/fqGSyGEEEII0ftoHcTjOYjLtROXeydu1y5c7p34/U7GjH6atLSbzn8Qg0jxLYQQQggheqxgsJ3W1j24XDtwuXficu3E7d5NKOQFQKlo7PbhDEidhT1+FImJ4w1OfG5SfAshhBBCiB4hEHDjcu/C5doeLrZdO/B4PkfrIAAWSzx2+2gyMu4k3j4ae/xo4mIvw2SKMjj5hZPiWwghhBBCdLtgsB2XezvNzVtoafkMl2sHbW2HOvZHRw8kPn40A1JnER8/hvj40cTEZKKUMi50J5DiWwghhBBCdKlQKEBb2xFcru00t2yhpbkCl3sXWvsBiLGmE58wlsFpt0QK7TFYrQMNTt01pPg2yJIlS1i5ciVmsxmTycSyZcv4wQ9+QG1tbceS8bm5uaxevZr77ruPAQMG8OMf/7jjvTU1NfzmN78x8hKEEEIIIU6hdZC2tiO4W/fS6t5Ha+s+Wj37aW09gNY+AEwmGwkJl5OV9Q0SE8aTkDAeq3WAwcm7jxTfBigrK2PNmjWUl5djtVqpr6/H5wv/D/niiy9SUHDqFJFPPPEE48ePZ8GCBSileO6559iyZYsR0YUQQgghAPD7G3G79+B27/77c+vejhshAWJiMoiLG47DMQN73HDs9lHExeVhMvXfErT/XrmBamtrSU1NxWq1ApCamnrO9gkJCSxZsoTvfOc7ADz++OMkJSV1eU4hhBBCCK1DtLUdxuXeFZnSbxdu9y7a2+s62kRFObDbR5KRcSf2uBHY7SOIjb0MiyXOwOQ9U/9e4fKtB6FuW+eeNC0f5j55ziZut5vp06fj8XiYNWsWxcXFzJw5k6KiolOGnXzlK1/hqaee6njf1KlTMZvNbNiw4ZzHlxUuhRBCCPFl/H1av52RYnsH7tY9BIMeAJQyExt7WWSmkZHY7aOwx40gOjq1198IeakudIVL6fk2gN1uZ/PmzZSUlLBu3TqKi4t58slwwX6mYScAVVVV1NXVoZTC7XZjt9u7O7YQQggh+pBAwIXLtQuXe8cZp/Uzm+3E20cxePCtxNvHYI8fSVzscMxmq8HJe7f+XXyfp4e6K5nNZoqKiigqKiI/P58VK1acs/3999/Po48+yq5du3jsscdO6REXQgghhDiXUMiHy7WT5uZymlsqzjCt34COaf3s8aOJt4/GZstEKZNxofuo/l18G2TPnj2YTCaGDx8OQEVFBdnZ2Wzfvv2M7d966y2OHTvGXXfdhcfjYdy4cSxatIjRo0d3Z2whhBBC9AJaB/F6a3C7d9HUXE5zczku1zZCofDkDv1pWr+eSIpvA7jdbhYvXkxTUxMWi4Xc3FyWL1/OrbfeyoIFCzrGfKemprJmzRq++93vsnr1apRSxMXF8bOf/YzvfOc7rF271uArEUIIIYRRgkEvLa5ttLr34mk7RJvnUPi5rbJj/mylokmIH8OQjH8iMXECiYlXYLUOMjh5/9a/b7jso/rDNQohhBD9jddbQ3NzOU3N5bQ0b8Hl3onWAQBMphhibdnYYnM6nuPicom3j5Ux2t1EbrgUQgghhOil2n31uFq20eLajsu1HVfLNtp9R4FwoZ2QMI6srG+SmDiBePtorNZBMj67l5DiWwghhBDCIFprvN5q3O6dkZlHduJybT9pDm1FbOwwkpOnkJAwjsTECdjtIzGZogzNLb48Kb6FEEIIIbqB1kFaPQdwtWyPTO+3E7d7F4FAS6SFidjYoSQnTSY+fmzkMRqLRaYX7kuk+BZCCCGE6GRah/B4DtLSsrVj6IjbvatjsRqTKQa7fSSDBt2A3T4qvGiNfQRms83g5KKrSfEthBBCCHGJfL4GWlo+o7llCy3NW2lxfdbRo20y2YiPDy9WkxDp0Y6NvQyTScqw/kj+qwshhBBCXIDW1gM0t5TT7q2jvb0Ob3tt+NlbRyDQFGllwm4fwcCB15GYMJ6EhMuJi8tFKbOh2UXPIcW3AcxmM/n5+WitMZvNLF26lMLCQg4dOsTQoUN5+OGH+c///E8A6uvrGTx4MPfccw9Lly7l0UcfxW6388ADDxh8FUIIIUTfFwi4OHp0DTW1r9LSsqVje1SUgxjrYGJiMkhMnIjNlklC/OXEx4/FYokzMLHo6aT4NoDNZqOiogKAd955h4ceeogPP/wQgGHDhrFmzZqO4vuVV15hzJgxhmUVQggh+hutQzQ2llFb+yrHjr9DKOQlLm44ubkPMiB1FlZrusydLb40Kb4N1tLSQnJycsfXNpuNUaNGsWnTJgoKCli1ahW33347NTU1BqYUQggh+p729qO43Xtp81bS1naYtrbKyOMIwaAbiyWewYPnRcZqX45SyujIog/o18X3Tz/5Kbuduzv1mCMdI/nBpB+cs01bWxvjx4/H6/VSW1v7hWXi77jjDl5++WXS0tIwm82kp6dL8S2EEEJ0Aq2DNDR8SFX1izQ0fAiEV/o2maKJicnEZsskKWkiiQkTGDBgNmZzjLGBRZ/Tr4tvo5w87KSsrIy77rqL7du3d+yfM2cOP/7xjxk0aBDFxcVGxRRCCCH6DJ+vnpqa1VTXrMTrrSY6egA5Of+KwzEdmy0Ta/RAWSFSdIt+XXyfr4e6O0ydOpX6+nqOHz/esS06OpqJEyfyi1/8gh07dvDGG28YmFAIIYTondraqmlsKqOhYT3Hj/8Nrf0kJU2OjN3+iqwSKQzRr4vvnmD37t0Eg0FSUlLweDwd27///e8zc+ZMUlJSDEwnhBBC9B7e9joaGz/qeHi9lUB4ZpKMjDvJyJiPPW64wSlFfyfFtwFOjPkG0FqzYsUKzOZT5/8cM2aMzHIihBBCnIXWGq+3iqamT2hs+pSmpk9oazsMgMWSQHLSZLIyF5KcPJW4uOEypET0GEprbXSGLlNQUKA3bdp0yrZdu3YxatQogxJ1j/5wjUIIIfofr7cWp3MjjY2lNDZ9Qnt7LQAWSyJJSVeSlHQlyclTiLePkkVtRLdTSm3WWhecr530fAshhBCiRwoG22hq+oQG5waczhJaW/cBEB2dSlLSJJKS7iE5aZL0bIteRYpvIYQQQvQYJxa4qa5+ifqG9wmFfJhM0SQlXsngwV/D4ZiBPW6EzLktei0pvoUQQghhOJ/PSW3dn6iufom2tkNERSWTnn4HqSlFJCVNwmy2GR1RiE4hxbcQQgghDBEMemlq+pS6utc5dvxNQiEfiYkFDBt6HwMGzJEl3EWfJMW3EEIIIbqF1prW1r00OEtwOjfQ1PQJoVA7ZrOd9PRiMtLvxG4fYXRMIbqUFN9CCCGE6DI+Xz1OZylOZwkNzg34fMcAiI3NJSNjPg7HdJKTJsuwEtFvSPFtALPZTH5+PlprzGYzS5cupbCwkEOHDjFq1ChGjBiBz+fjqquu4n/+538wmeQObiGEEL1DKNROU3M5zsgMJS7XDgAsliQcjkJSHDNwOKYTE5NucFIhjGFo8a2USgKeA8YCGvg6sAdYBeQAh4DbtdaNKnxb86+A6wAPsFBrXW5A7Etms9moqKgA4J133uGhhx7iww8/BOCyyy6joqKCQCDANddcw+uvv868efOMjCuEEEKcldYaj+cgTud6GpwbaGr6mGDQg1IWEhMnMGzY90hxzCA+fozMvS0Exvd8/wp4W2t9q1IqGogFfgi8r7V+Uin1IPAg8ANgLjA88pgM/Dby3Ku1tLSQnJz8he0Wi4XCwkL2799vQCohhBDi7Pz+ZpyN4aEkzoYSvO01ANhsOaSlzSPFMYPk5MlYLPEGJxWi5zGs+FZKJQBXAQsBtNY+wKeUuhkoijRbAXxAuPi+GfhfHV6S8yOlVJJSarDWuvbLZqj7r/+ifdfuL30NZ2IdNZK0H/7wnG1OLC/v9Xqpra1l7dq1X2jj8Xh4//33efzxxzs1nxBCCHGxQqEALS0VOJ0baHBuoKVlKxDCbLbjcBSS7biXFMcMbLZMo6MK0eMZ2fM9DDgOPK+UGgdsBu4HBp0oqLXWtUqpgZH2GUDlSe+vimw7pfhWSt0N3A2QlZXVpRfwZZ087KSsrIy77rqL7du3A/D5558zfvx4lFLcfPPNzJ0718ioQggh+qm2tiMdK0s6naUEg27ARELCOIbm/CsOx3QSEsZjMhn9R3Qhehcj/8VYgAnAYq31x0qpXxEeYnI2Z1rKSn9hg9bLgeUABQUFX9h/svP1UHeHqVOnUl9fz/Hjx4G/j/kWQgghulMg4Kax8aPINIAltLUdBsBqHcyggdfhSLkKR/JUoqKSDE4qRO9mZPFdBVRprT+OfL2acPF99MRwEqXUYODYSe1P/nvWEKCm29J2kd27dxMMBklJScHj8RgdRwghRD8RCgVwuT6jwbkRp3MDLS0VaB3AZLKRnDyFzCF34XBcRWzsUFnKXYhOZFjxrbWuU0pVKqVGaK33ANcCOyOPfwaejDz/OfKWvwDfUUq9TPhGy+ZLGe9tpBNjviF8l/iKFSswm+UOcCGEEF1Ha01b2yEanBtodG6ksekjAgEXoIiPH0NW5jdwpMwgKXECJpOsLClEVzF6oNZi4MXITCcHgEWACfijUuobwBHgtkjbNwlPM7if8FSDi7o/bucIBoNn3J6Tk9Mx9lsIIYS4VH5/C42NZZEFbkrweqsAiIkZwsCB10UWuJlCdLTD4KRC9B+GFt9a6wqg4Ay7rj1DWw38a5eHEkIIIXqZUChAe3sNnrYjtLUdoc1ziOaWLbS0bEXrIGazneTkKWRn3Y3DMQ2bLVuGkghhEKN7voUQQghxkbTWNDaWUl2zCpdrO15vNVoHOvabTNHY40aSnf0vOBwzSEwYj8kUZWBiIcQJUnwLIYQQvYTf30xt3Z+orn4Rj+cgUVHJJCdPZeDA64i1ZWGzZWGzZWO1DkIpk9FxhRBnIMW3EEII0YNpHaKl5TOqa17m6NE3CIW8JCZcwejRv2DggLmYzXJzpBC9iRTfQgghRA/T3n4scpPkBpzODfj9TkwmG2lpX2VIxgLi40cbHVEI8SVJ8S2EEEIYLBj00tS8CWfDepzODbhb9wAQFZVCiuMqHI7pDBgwC4sl3uCkQohLJcW3QZYsWcLKlSsxm82YTCaWLVvGp59+yi9/+Us+//xzjh8/TmpqqtExhRBCdAGtNa2te3E6N9DgLKGp6RNCoXaUiiYpqYDctP/A4ZiB3T5Sxm4L0cdI8W2AsrIy1qxZQ3l5OVarlfr6enw+H9HR0dxwww0UFRUZHVEIIUQn8/mcOCPDSJzODbT7jgIQG5tLRsb8yJzbkzGbbQYnFUJ0JSm+DVBbW0tqaipWa/gmmRM93Onp6UbGEkII0YlCIR/NzVs6FrhxuXYAGoslEYdjGimOGTgc04mJkZ/9QvQn/br4LvnjXuor3Z16zNRMOzNuzztnm9mzZ/P444+Tl5fHrFmzKC4uZubMmZ2aQwghRPf6+/LtJTidG2hs/IhgsBWlzCQkXMGwod/FkTKDhPixKGU2Oq4QwiD9uvg2it1uZ/PmzZSUlLBu3TqKi4t58sknWbhwodHRhBBCXISzLd9ui8kiLe2rpDimk5w8VW6UFEJ06NfF9/l6qLuS2WymqKiIoqIi8vPzWbFihRTfQgjRw4VCAVyuzyJTAJacZfn26cTGZhsdVQjRQ/Xr4tsoe/bswWQyMXz4cAAqKirIzpYf1EII0VMEg17avJW0tR2hzXOYtrYjeNoO0dKylUCgBVAkJFwuy7cLIS6aFN8GcLvdLF68mKamJiwWC7m5uSxfvpxnnnmGn/3sZ9TV1XH55Zdz3XXX8dxzzxkdVwghup0/6Mft79x7ck4XCLTQ7q2i3VuFL/J84uH3HTulrckchzVmCImOIhKSphKfeCWWqKSO/c2+rs0qer4oUxRxUXEopYyOIno4pbU2OkOXKSgo0Js2bTpl265duxg1apRBibpHf7hGIUTv0dzeTJW7Cl/Q94V9WmucXieVrkqOuI5Q6aqkylVFbWstIR26pPMqNAlmTYpFk2rRpJpD4ecoTYo5RNxp9zw2B6EhYKI+oKiPPDdEXreGwkcU4lyiTdGk2FJIiUkh1ZZKii2FRGsi5l50g63NYvvCI8YSc9ZrMJvMRJuiiTJHhZ9NUUSZo87aXimFRVmwmP7+MCtzn/jQopTarLUuOF876fkWQoh+yhvwUuWqotJVSWN7Y6ccMxAKUNdaR6WrsqOgdvlcF/TeJGsSWfFZjBswjhuG3UByTDLqRMGr/ZgDjZiDjZhCZ+5lVtqPOdiIOeDEEnRiDjhRBDr2a0wEzYkEzQ6ClhRcZgdBSzJBs4OA2QGmaOKAOEAGAoqL5Qv6cHqdNHgbqG+rp7a1lu0N22lqb4Je0s8ZInTJH3q/LIuyYFImTMqEUqqjIDcrM6bIQlMmZcJEeL9JmVCojtcnvjYpE/dPuJ9rsq4x5DouhBTfQgjRS4R0iGOeYx29w2frTT4brTXNvuZwYdxSybG2Y+d/05dgVmbS7elkxmdyXep1ZMZnMsQ+BFvUmRePSbImkRmfSXx0fGS6vsM4nRtocW0Lj7luO0x7+1EupIIxmWKw2TKxxV9OrC0bmy2r4xETkyHjsoU4B601/pCftkAbbYE2PAFP+LW/DX2Wf3+BUAB/yI8/6Mcf8uML+fAFfWct4rXWBHSg432BUPh1UAcJ6RBa647XIR0iqINA+OefRqO1PmW71rrjQ8OJfQnRCV3zDeokUnwLIcRpTi5yTzyqXFV4A15D8vi1nzp3HVXuKtqD7R3bzcpMtDn6oo4VFxVHZnwmU9KnkBmfSWZ8JlnxWaTaUjvlz74KhcPmIOoiitxAwIWzsZRq5wYaGkrweisBiIpKITY2m+TkKdhsOcRGiujo6IFnzKpM0URHpfSJP18LYQSlFNHmaKLN0SRaE42O02dJ8S2E6LWqXFWU1pRSVlPG/qb9Z+2ZuRjBUJBjnmP4Qn/vUbYoC4Ptg7FH2S/5+F+GUorshGymZ0wPF8wJ4aJ5cNxgLKbe8WNc6xDt7XV42g5HerOPdMwk4m7dHZmuL47k5KlkZ30Th2M6Nlu2FNJCiD6nd/zUFkIIoNXfyqd1n7KxeiNltWUcbjkMQFpcGvmp+Z1SiCoUA2MHhodKxA/pdUVuT+LxHMLp3ECDsySy2uPfx2orFUVMTAaxsdlkp9yDI+Uqma5PCNEvyG8TIUSPFdIhdjXsorSmlI01G9l6bCsBHcBmsXFl2pXcOfJOCtMLyUnIkR7SHiAQcNHYWBZeXr1hA23eIwDExAxh0KAbSIgfGxl/nU1MzGBZYl0I0S9J8W2QJUuWsHLlSsxmMyaTiWXLlvHMM8+wadMmoqKimDRpEsuWLSMqSnqBRP9ytPUoZbVllFaXUlZbFp4pABjlGMVdY+6iML2QKwZecdFjnUXn0zpIS8u2yNLqG2hp2XLK8JHMrK+T4piOzSanTdEOAAAgAElEQVQfjoQQ4gQpvg1QVlbGmjVrKC8vx2q1Ul9fj8/nY8GCBfzf//0fAPPnz+e5557j3nvvNTitEF2rLdBG+dFyNtZs7Bi7DZBqS+WqIVcxNX0qUwdPJcWWYnBSAeD11nQMJXE6NxIINAOK+PixkaXVZ5CYeAUmk3w4EkKIM5Hi2wC1tbWkpqZitVoBSE1NBSA9Pb2jzaRJk6iqqjIknxBdLRAKsL5qPav3rubj2o/xhXxEm6KZMGgCN112E4XpheQl50lvaQ8QDHpobPw4UmxvwOP5HABr9CAGDPgKDsd0HMnTiI52GJxUCCF6h35dfK/7w3KOHT7QqcccmD2Mqxfefc42s2fP5vHHHycvL49Zs2ZRXFzMzJkzO/b7/X5eeOEFfvWrX3VqNiGMVt9Wz6t7X+WVva9w1HOUgbEDuX3E7UzLmMbEQROxWc48D7ToWloHO2Yf8UTm1e742nMIrf2YTFaSkiaRkV6MwzGDuLjh8uFICCG+hH5dfBvFbrezefNmSkpKWLduHcXFxTz55JMsXLgQgG9/+9tcddVVzJgxw9igQnQCf8hP+dFyXtn7Cu8ffp+ADjB18FQemvQQMzNnyiwiBvF6ayPDR9bjdJYSCDR17AsvVBO+MTI15RqSHYUkJV6J2Ww1MLEQQpyf1uEpZ3ty50C//q13vh7qrmQ2mykqKqKoqIj8/HxWrFjBwoULeeyxxzh+/DjLli0zLJsQl6rSVUlpdSmlNaV8UvcJbr+bhOgE5o+az215t5GTmGN0xH7F72/p6M1ubtlCQ0MJHk94bH109EAGpF5DUtIkbLE5xNqyiY4e0KN/cQkhxOl8bR52bfiAre++xTWL7mHIqLFGRzqrfl18G2XPnj2YTCaGDx8OQEVFBdnZ2Tz33HO88847vP/++5hMJoNTCnHh3D43n9R9QmlNuOCudIVXKEyPS2fO0DkUphcyI2MGMZYYg5P2bVoHcbl24HRuwN26NzJs5PBpvdrh4SPp6beR4phBXJyMrRdC9F5HD+xn63tvsXvDh/jbvQzIHkowEDA61jlJ8W0At9vN4sWLaWpqwmKxkJuby/Lly0lLSyM7O5upU6cCMG/ePB555BGD0wrxRcFQkJ0NOzuK7a3HtxLUQWwWG5PSJrFg1AKmpU8jO0FWKOxq4eEjG2lwrqexsRS/vxGAmJhMYm1ZDBo4NzyEJDYbmy2bWFsOZrN8CBJC9F5+r5fdpevZ+u5bHD2wD0u0lRGFMxg3ay5puT2/Q0GKbwNMnDiR0tLSL2wP9PBPaqJ/q2ut6yi2P6r9iOb2ZgBGp4xm0dhFFKYXMn7AeKLMMjd9VwsEXNTWvU5NzSrc7l1AePhIasrVOBwzcDgKiY5ONTilEEJ0roaqI1T87U12rl+Lr81DypAsrll0D6NmXE1MnN3oeBdMim8hxBm1BdrYVLepo+A+0ByeGWigbSAzh8xkWvo0pqRPwREjU8x1F5drF9XVL1J39M8Egx7i4/PJzX2QFMdVMnxECNEnBQN+9n1cytZ336Jq13bMFgvDJ09j3OzryBgxulf+3JPiWwgBhO8Q39u4l401GymtKaX8aDn+kB+r2UrBoALmDZ9HYXohuUm5vfKHXW8UCLhpazuCy72TmppVNDeXYzJZGTToRoZkLCAh4XKjIwohRJdorK1m27p32fHBe3iam0gcOIgZ8xcy9uqvEJuQaHS8SyLFtxD9WH1bPWU1ZZTWlFJWU0aDtwGA4cnDmT9yPoXphUwYNEFulDxNKOTD663GE5lBJBBwdd5x26oixz2M3+/s2Gez5TA890cMHjyPqKikTjmfEEL0JH5fO/s/LuWzte9QtXM7ymRi2IQrGTdrLjnjJqD6yGQUUnwL0Y+0B9vZcmxLeChJdSl7GvcAkGxNZkr6FKalT2Nq+lQGxg40OKnxTvQ6t0UWnfGctPCM11sLhLrgrIoY62BstiwGpM7CZsvGFptFrC0Hu30kSvWNXzxCCHGy+iOH+Oz9d9hZspb21lYSB6Ux/Y67GDPzWuyOFKPjdTopvoXow7TWHGg+0DFue1PdJrxBLxaThSsGXsH9E+6nML2QkY6RmPpoYad1iPb2o18opL3eWtDBL7YniNdbh9/fcMr2qCgHNls2SYkF2NKyIovQhBeiiYpKBC59KI5SJpQyX/JxhBCipwsGAny+6SMq3vkrlTu3YY6KYvikQvKv+QcyR4/tM73cZyLFtxB9THN7M2W1ZR2L3Bz1HAUgJyGHW4bfwrT0aRSkFRAXFWdw0s5z+jCQ0x+hUHtHW6XMxMRkEBMzBJPpTDOzKOLtY7DF5mCzZREbKbItlvjuuyAhhOijWpsa2fb+O2x97y3czgYSBgziqgWLGHv1V7DFJxgdr1tI8W2QJUuWsHLlSsxmMyaTiWXLlrF8+XI2bdqE1pq8vDz+8Ic/YLf3nqlzhHG01mw5toWX97zMu4ffJRAKEB8dz5TBUyhML2Rq+lQy7BlGx7xkbW3VuFzbO1ZrPFFse701nDwMxGSyEWsLD9dIcVwVGb6RTawtC6s1HZMsaS+EEN1Gh0JU7tzGZ++/w76PSwkFA2RffgWzvvlthl5RgMnUv/7iJ7+BDFBWVsaaNWsoLy/HarVSX1+Pz+fj6aefJiEh/Knve9/7HkuXLuXBBx80OK3oyVr9raz5fA0v73mZ/U37iY+Kp3hEMXOHzmVMyhgsvbzIDARaaWr6mAbnepzODXg8Bzv2RUUlY7NlkZh4BWlpN0d6qcO91bI8uhBCGM/d6GTHB++xfd27NB2txRoXx7jZcxk/+3oc6UOMjmeY3v2buZeqra0lNTUVq9UKQGrqqYthaK1pa2uT4kGc1d7Gvfxxzx954/M38AQ8jHKM4rHCx5iTM4fYqFij431pWoc6lkdvcJbQ3FyO1n5MJhvJyZPJyFhAUmIBsbE5MgxECCF6IB0KcXDrZj57720OlH+KDoXIHJ1P4W3zyZ1cSFS01eiIhuvXxXfTG5/jq2nt1GNGp8eRdONl52wze/ZsHn/8cfLy8pg1axbFxcXMnDkTgEWLFvHmm28yevRofvGLX3RqNtG7+YI+3jv8Hqv2rKL8WDnRpmjmDJ1D8Yhi8lPze+2Htfb2ozQ4S3A6N+B0buyYXs9uH0VW5iIcjukkJRVgMskPbCGE6Kn8Xi871q+l/M0/01hbTWxiEgU3ziP/6q+QPLj3D3vsTP26+DaK3W5n8+bNlJSUsG7dOoqLi3nyySdZuHAhzz//PMFgkMWLF7Nq1SoWLVpkdFxhsGp3Nav3ruZP+/6E0+skMz6TBwoe4ObLbiYppvfN9xwMemlq+hSns4QGZwmtrXsBiI5OJSXlqvDy6MnTsFoHGJxUCCHE+bic9VS8vYbP3nsbb6ubQcOGc93iB8ibMh2zRcrMM+nX35Xz9VB3JbPZTFFREUVFReTn57NixQoWLlzYsa+4uJinnnpKiu9+KqRDbKzeyKo9q1hftR6lFDOHzOSOEXcwJX1Kr5oWUGtNa+vejt7tpqZPCIXaUSqapKQCBqfdgsMxA7t9hMxjLYQQvUD4BsrtbFv7Dns/2oAOaXKvnMKE62/utUu+d6d+XXwbZc+ePZhMJoYPHw5ARUUFWVlZ7N+/n9zcXLTWvPHGG4wcOdLgpKK7NXobeW3/a/xxzx+pdleTEpPCN/O/yW15tzHYPtjoeBfM52vA6dwY6d3egM93DIC4uOFkZMwnxTGDpKRJmM02g5MKIYS4UC3Hj7Hjw/fZ/sF7tBw/SrQtlvGzr+eKuTeRNCjN6Hi9huHFtwqvKLEJqNZa36CUGgq8DDiAcuCftNY+pZQV+F9gItAAFGutDxkU+5K43W4WL15MU1MTFouF3Nxcnn32WW655RZaWlrQWjNu3Dh++9vfGh1VdJP9jfv5/fbf886hd/CH/BQMKuC7E7/LtZnXEmU+01zUPYPWGp+vvmPqv9bW/TgbN+By7QDAYknC4ZhGimMGDsd0YmJ6zwcIIYQQ4HW7ObjlU7Z/+D5Htm8FIGvM5Uy/45/InTRVbqD8EgwvvoH7gV3AiZnVfwo8rbV+WSn1LPAN4LeR50atda5S6o5Iu2IjAl+qiRMnUlpa+oXtGzduNCCNMNIxzzF+U/EbXt//OjaLjVvzbuX2vNvJTc41OtoXnBg+4nRuoKl5c0fBHQx6OtooZSEx4QqGDf03UlKuIj5+jKzYKIQQvYgOhTh26AAHKzZzsGIztft2o0MhEgcOovDW+YyZeS0JAwYaHbNXM7T4VkoNAa4HlgDfU+FBQtcA8yNNVgCPEi6+b468BlgNLFVKKa217s7MQnSGVn8rz29/nv/d+b/4Q34WjFrA3fl397gbKH0+J87GjTgbwuO1233h1TJttmxiY4eRnDTllGXWbbYMmZVECCF6oaMH9rPl7TUcrNiEp7kJgEHDcpn81dvIGV9A+vARfXrJ9+5kdM/3L4H/AE5M2JsCNGmtA5Gvq4AT89NkAJUAWuuAUqo50r6+++IKcWn8IT+v7XuN31T8BqfXyZycOdw34T4y4zONjgaEl2lvbt7SMVbb5doOaCyWxNOGj6QbHVUIIUQnOHboAKWvrOTzTR8RbYtl6BUFDB0/kZxxE4hLSjY6Xp9kWPGtlLoBOKa13qyUKjqx+QxN9QXsO/m4dwN3A2RlZXVCUiEu3THPMV7d9yqr967mmOcYEwZOYOk1S8kfkG9orvCCTodocG7A6SyhsfEjgsFWlDKTkHAFw4bejyPlKhLix8rwESGE6EOOHzlE2Ssr2fdJKda4OApvX8CEuTdjje29C7X1Fkb2fE8DblJKXQfEEB7z/UsgSSllifR+DwFqIu2rgEygSillARIB5+kH1VovB5YDFBQUyJAUYRitNZ/WfcrLe15m3ZF1BHSAaenT+MnUnzAjY4ZhUzH5/S00NpZ1zLPt9VYBYIvJIi3tq6Q4ppOcPFVWkBRCiD7o6IH9fPKXV9lbVkK0LZapt97JhOtuJibObnS0fsOw4ltr/RDwEECk5/sBrfUCpdQrwK2EZzz5Z+DPkbf8JfJ1WWT/WhnvLXqiQCjAa/tf44WdL3Cw+SCJ1kT+cfQ/clvebWQldP9fY0KhAC7Xtkjv9npaWraidRCz2U5y8hSys76FwzGd2Nicbs8mhBCi6/m9XnaXrmfru29x9MA+omJsTJlXzITrv4rNLh0t3c3oMd9n8gPgZaXUE8AW4PeR7b8HXlBK7Sfc432HQfmEOCOtNesq1/H05qc51HKIsSljeWLaE/xDzj8QY4np1ixtbdUdPduNjaUEAi2AIiE+n+zsf8HhmEFiwnhMpp47jaEQQohLU3/kEFvfe5ud69fia/OQmpnNNV//F0bPuBprbJzR8fqtHlF8a60/AD6IvD4ATDpDGy9wW7cG60JLlixh5cqVmM1mTCYTy5YtY/LkyQAsXryY559/HrfbbXBKcaG2Hd/GLzb/gs1HN5OTkMOvrv4VV2de3a1DS/z+FmrrXqWmZhWtrfsAsFrTGDhgDg7HdByOQqKi5OYZIYToy/y+dvZ9tJGt771NzZ6dmC0W8qbOYNysuaSPGCWrT/YAPaL47m/KyspYs2YN5eXlWK1W6uvr8fl8AGzatImmpiaDE4oLVdlSyTNbnuHtQ2/jiHHw8OSHmZc3j6hu7FFucW2nuupF6o7+hVDIS0LCFQwf/jApjhnExl4mP2iFEKIfaKiq5LP33mLn+rV4W90kD07nqn/8OmNmXktsQqLR8cRJpPg2QG1tLampqVit4fmQU1NTAQgGg/z7v/87K1eu5LXXXjMyojiHkA6xsXojq/asYn3VemIsMdxz+T0sGruIuKju+TNeMOjl2LG/UlW9kpaWCkwmG2lpNzMkYwHx8WO6JYMQQghj+do87PukjG1r/0b17h2YzBaGT5rK5bPmkjkmXzpfeqh+XXy/9dZb1NXVdeox09LSmDt37jnbzJ49m8cff5y8vDxmzZpFcXExM2fOZOnSpdx0000MHixLcPdEjd5GXtv/Gq/seYUqdxUpMSl8M/+b3DHyDgbGds9qXx7PIaqrV1JT+yqBQBOxsZeRN/wR0tJuISoq4fwHEEII0asFAwEOf7aFnSXr+HzTxwR87SQNGsyM+QsZWzSL2MSetVib+KJ+XXwbxW63s3nzZkpKSli3bh3FxcXcd999vPnmm3zwwQdGxxOncXqdPL35ad488Ca+kI+CQQXcP+F+rs26lihz1w8vCYUCNDSso6r6RZzOEpSyMGDAbDIy5pOcNEV6NoQQoh84emA/2z94lz2lJbS5WoiJT2DMzGsZNeNq0vNGyu+CXqRfF9/n66HuSmazmaKiIoqKisjPz+fOO+8kJSWF3NxcADweD7m5uezfv9+wjALWHVnHo2WP4vK5+Nrwr1E8opjc5NxuOXcoFKC2djUHDy2lvb0WqzWNYUO/S3p6MVZr9/S0CyGEMFZDVSUbXl7B/k8/whIVzWUFkxk1o4iccRMwW2TGqt6oXxffRtmzZw8mk4nhw4cDUFFRwT333MPSpUs72tjtdim8DeT2ufnZpz/jtf2vMSJ5BL+b/TvykvO65dxaaxoa1rH/85/R2rqPxMQJ5OX9mNSUazGZ5J+sEEL0B66GekpfWcmOD94jKsYaWYHyJpkisA+Q3+QGcLvdLF68mKamJiwWC7m5uSxfvtzoWCLi07pPeXjDw9R56vhW/re4d9y93TK8BKCl5TP27X+SpqaPsdlyyM//HwakzpY/JwohRD/hdbv55M+vsOWtNwiFQlwx90Ym33K7zFjSh0jxbYCJEydSWlp6zjYyx3f3aw+280z5M7yw8wUy4zNZMWcF4weO75Zzu1v3cejQbzh69A2iohzk5T1KRvodsgiOEEL0AwG/n8OfbWFPWQmfb/oIn9fL6OlFFN7+jyQOHGR0PNHJpPgWAtjRsIMflfyIz5s/p3hEMd+b+D1io2K79JyhkI/jx/9GVfVKmpo+xmSykpN9L9nZ92CxyHK/QgjRlwUDfg5vq2Bv2Qb2f/oR7Z5WYuLs5E2ZzhVzbmRgzjCjI4ouIsW36NcCoQDPbXuOZVuX4Yhx8OysZ5mWMa1Lz+n11lBd8zI1Navw+eqJickk97L/YPDgW4mOTunScwshhDBOKBSkaud2dm/8kH0fl+JtdWONjSP3yinkTZ1Odv54uYmyH5DiW/RbB5sP8qMNP2Jb/TauG3odP5z8QxKtXTemrq2tigMH/pu6o28AmtSUq8kYsoAUx1UoZeqy8wohhDCO1prafXvYXfohe8s20NrUSJQ1hssKJjNy2lVkXz4BS5QU3P2JFN+i3wnpEC/tfomnNz9NjCWGp2Y+xZycOV12Pr+/mUOHfkNl1QsopcjK+jpDMv4Jm21Il51TCCGEsULBINvW/o1P/7Ka5mNHMUdFMXR8ASOnzWTYhAKirDFGRxQGkeJb9Csev4eHNz7Mu4ffZUbGDB4rfIwBsQO65FyhUDtVVf/HwUO/IRBoYXDaPIYN+zdiYmQFUyGE6MsOVmzmwxd+T0PVEdLzRjH11vnkXjlFpgkUgBTfoh+pdddy37r72Nu4lwcKHuCu0Xd1yRR+wWA7R4+9wcGDS/F6K3E4ZpCb+yDx9pGdfi4hhBDGCnm9BJ1Ogm43DXU1bHznDSr37yHRkcqcOxYybHQ+Sil0VTVeo8OeTyhEyOMh1NpKyO0OP7e2EvJ40Fpf+HG0hmAI7fejA4Hww+9DBwIQOstxQiF0KBh+38nPIQ2hEGiN1iHQRL4OhTNF9p+8L/Xb92KfMaMzviNdQopvgyxZsoSVK1diNpsxmUwsW7aM3/72t3z44YckJobHHf/hD39g/Pjumequr9tybAvfXfddfEEfS69Zyowhnf+P0uM5THXNSmpqVhMINGG3j2L8+BWkOKZ3+rmEEOJChDwefJWV+I4cwX8k/BxsbDQ6Vq8W8rYRdDYSdDoJNDai29pojbZwYGAylY54ooIhRh1tJPuzzzGt+5hDRgc2UlQUymJBnXi2WMBsPnNbBcpkBrPp1GeTCUwqfG+UUmAyhTvOIq9P7DvRVpkU6mzn6CGk+DZAWVkZa9asoby8HKvVSn19PT6fD4CnnnqKW2+91eCEfctr+17j8Y8eJz0unV/P+TXDEjtv+qZQKEBDwzqqql/E6SxBKTMDUmeTkTGf5OSpsjiOED2YDoUIHD+O/8gRfEcq8VVGCtSqKnRbm9HxLpEm0NRE8Hj9KVvNiYmYB6TKz6ZLoKKtmB0O/FlDqNF+DrsacbqaMZlMjB15OVeMv5KY3jqe26QwxcZhiovFFBeH2W7HFBeHyWY7e9F81mOZ5P+zs5Di2wC1tbWkpqZitVoBSE1NNThR3xQIBfjvzf/NCztfYMrgKfx85s87bTaT9vbj1NSsorrmZdrba7FGD2Lo0O+SkX47VqssiCBET6H9fvzV1af2/lZW4jtyGH9lFbq9/e+NzWai0tOJzszENKj3/zuOSYgnOjOL6OwsojKziM7KxJyQYHSsXs3tbGDXxg/ZU1rC0QP7ABicO4KZX72VvCnTSUjtmnuIRN+iLmoMTy9TUFCgN23adMq2Xbt2MWrUKAD27v1PXO5dnXrOePso8vJ+fM42breb6dOn4/F4mDVrFsXFxcycOZOFCxdSVlaG1Wrl2muv5cknn+wo0C/GydfYX1W6Knl4w8OUHytn/sj5/PuV/47FdGmfNbXWNDV9QlX1/3H8+N/QOoAjeRoZGQtITb0W0yUeXwjx5QUaGmjbsgXf4cP4jlTirwz3ZvtrasLjQyNUTAzRmZlEZWURnZl5SmEaNXgwSqZ8E2fQUn+MT15fzfZ1fyMYCDBo2HBGTJ1O3pTpsgKl6KCU2qy1LjhfO6kWDGC329m8eTMlJSWsW7eO4uJinnzySf7f//t/pKWl4fP5uPvuu/npT3/KI488YnTcXkVrzep9q3nq06cwKzP/Nf2/uPGyGy/xmCFqalZRWbWC1tZ9WCwJDBlyF0My5hMbO7STkgshLkbI56OtvJzWjRtxb9xI+86/d6SYk5KIysrCNm4cCTfeEO79zQoX3JYBA+RP4eKCNR2t45PX/8iOD9cCMKboWq68cR7JgzMMTiZ6s35dfJ+vh7ormc1mioqKKCoqIj8/nxUrVrBw4UIArFYrixYt4uc//7lh+Xqj457jPFL6CBuqNzB58GSemPYEaXFpl3TMYNDDzp3/wbHjbxEfP5ZRI3/KoEHXYzbbOim1EOJctNbhcdmVlfgOH8FXeQTvzp14Pvk0PC7bYiH2iisY8G//RtyUyUQPG4Y5Pt7o2KKXa6yr4eM/rWJnyTpMJhP51/4Dk27+GgmpA42OJvqAfl18G2XPnj2YTCaGDx8OQEVFBdnZ2dTW1jJ48GC01rz++uuMHTvW4KS9x9uH3uaJj56gPdDOg5Me5M6Rd2K6xFUjvd4atn52D273LnJzHyIr8xvSYyZEF9B+P/7a2lOGi4THZx/54s2PJhPR2dkkzZtH3LRpxE6ahNkucyeLzhEKBdn81z+zcdULKBRXzLmRK2+ch92RYnQ00YdI8W0At9vN4sWLaWpqwmKxkJuby/Lly7n99ts5fvw4WmvGjx/Ps88+a3TUHq+5vZklHy/hrYNvkZ+az5LpSxiaeOlDQZqaN/PZZ/cSCrUz7vLfkZp6dSekFaL/CrW14aus/OLMIkeOhMdlB4MdbZXVSlTmEKIzs4grLCQqK5PoyBjtqPR0VHS0gVci+ipnTTVv//Zpavfu5rKCKcz6xr1SdIsuIcW3ASZOnEhpaekXtq9du9aANL3XxuqNPLLxEZxeJ98e/22+lf+tS76pEqC29lV27X6YmJg0xl2+kri43E5IK0TfprUm2NQUHh5yJDKbSGRmEf+RIwSOHz+lvSkhgeisLGz5Y0m47jqiIwV2x7hs06X95er/s3ff8W1WVwPHf4/28t4jtmNn7ziLDBJW2KVlBmhpC4X2ZRcCZUOBQilQKHs0UKAtSRgFyl6BkL23M+043vGQLGtY87nvHzImIUlJkBx53C8ffWRbjx8fE0s6ujr3HEk6XEJVWfvx+yye+yo6g4HTrpnN0GnHyXc6pS4jk2+px/EGvTy25jHmb59PSVIJT574JMPThkd9XlUNUl7+CFXVL5GSMoWRI55Cr0+OQcSS1LsEqqvxrlwV6Syyzwq26nLtd5wuMxNDQQHWY4+NbHjs990KtjZZ3rek+HM01PHpc09Qu20LxaUTmHnFNXK1W+pyMvmWepT1jeu5ffHt1Lhq+OWwX3Jd6XUYtUfejvH7PJ5yyspuos21kfy8Sxg48A40GtlyTJIAwm433hUrcC9ejGfJUoJVVZEbdDr0ebkY+hWQNHpUZ8s+Q0EB+vz8yGAOSeqG2l1trHjnDdZ/9iE6vYFTr7qBYdNPkKvd0lEhk2+pR1CFynMbnuPFjS+SbcnmpVNeYkL2hKjPK4RKTc1r7Cp/GK3WwogRT5GVeXoMIpaknkf1+QjW1BCoquocSOPbvp329eshHEaxWLBOmkTqL3+JdfIxGAoLI+OiJamHCPp8rPnoPVb9922CPh/Dpp/A1At/QUKqHHYnHT3yUVPq9jxBD7cuupWvq7/mrJKzuG3ibdgMtqjP6/PVUbb1Dzgcy0hLO56hQx7EaJRtpKTeLex07tNVpGMQTVUVgepqQnv37nesJiEBQ//+pP3mN1inTcUyZozc7Cj1SOFQiE0LPmP523PxtDooGX8M0y68hPR+hfEOTeqDZPItdWs1rhquXXAtu527uXXirVw85OKo3xYUQtDQ8A7bd9wLCIYMeZDcnAvk241SrxR2OvEsW45nyRI8S5ZEOovsQ5uRjqGgEOvkyZGuIvsMpNEmJ8v7hdTj7V6/hq9eeQFHfR15Q4bzkxtvJ29w354CLcWXTL6lbmtVwypu/PpGwiLMsyc9y5TcKVGfMxBoYdv2O2lq+ozkpAkMG/YwZnNBDKKVpCKrpqQAACAASURBVPgT4TChhgYCVVV4V6/Bs3gx7Zs2gaqisdmwTj6GlJ9fjKGwMFKf3S8fjcUS77AlqUu47M18/eocdixfTEpuPmffcg/9x46XLyiluJPJd5w88MADvP7662i1WjQaDS+88AITJ07kzjvv5M0330Sr1XLllVdy3XXXxTvUuHhzx5s8uPxB8hPyeeqEpyhKKor6nE1Nn7N12x2EQi4GlNxCQcFvUBRt9MFK0lGk+v2dddn7Tn0MVlUTrK1FBIORAzUaTCNHkP5//4d12lTMI0ei6OUmYqn3U8Nh1n/6AUve+BdqKMzUWZcw/ifnoJN//1I3IZPvOFi2bBkffPABa9euxWg00tzcTCAQ4JVXXqG6uppt27ah0WhobGyMd6hHnS/k49HVjzJ/+3ym5k3l4ekPk2hIjOqcoZCLHTv/RH39W9hswygd+09stsExiliSYi/scn034XHfgTTV1YQaGkCIzmM1Viv6ggKMgwaRMPOkznZ+piFDZDs/qc+p37WdL/7+LI2V5RSNGceJl/4fydk58Q5LkvYjk+84qK+vJz09HaMx0iIvPT2yy/q5557j9ddfR9MxXCIzs29t/tvSsoXbF91OhbOCXw37FTeMuwGtJrqVaYdjOWVb/4DPV09R4VX0738tGo3cMCZ1H6rfT/uaNbiXLKF99RoCe/YQbm3d7xhtWlqkX/bECZ3t/L5NsrWpqfJtdKnP83u9LJ73Kus/+whbcgo/ueFWBk6aKu8bUrfUp5Pvu3bWsNndHtNzjrCZuX9g/v885uSTT+a+++5j0KBBnHTSScyaNYsZM2ZQXl7O/Pnzeeedd8jIyODJJ59k4MCBMY2vOwqqQeZsmsOLG14k1ZzKCye9wJS86Oq7hVCp2P0ElZVPYzYXMX7cGyQljY1RxJL04wkhCJSX41myBPfiJXhXrUL4fKDXYx49ioSTT8ZQWNCZXOvz+6G1WeMdtiR1WztXLmXBP17A7bAz9pQzmTrrEoxyL4PUjfXp5DtebDYba9asYdGiRXz11VfMmjWLhx56CL/fj8lkYvXq1fznP//hsssuY9GiRfEOt0vtdu7m9kW3s7llM2cUn8FtE28jyZgU1TlDIQ9lW2+iqekzcnLOY/Cge9Bq5QOxFF+qx4Pz/Q9wzJ2Lf/t2AAz9+5N8/vlYp07BOmECGqtMsiXpcLlamlnwj+fZtWo5GQVFnDX7dnIGyJJCqfvr08n3D61QdyWtVstxxx3Hcccdx8iRI3n11VfJz8/n3HPPBeDss8/m0ksvjVt8XU0VKnO3zeXxNY9j0pl4dMajnFJ0StTnbW+vYeOm3+F272DQwLvIz/+VfNtRiiv/zp045s7D+d57qB4PxqFDybr7LhJmzECflxfv8CSpx1HVMBs+/5jFc19FDYU59uJfM+6Mn6GVA5+kHkL+pcbB9u3b0Wg0nSUl69evp7CwkOHDh7NgwQIuu+wyFi5cyKBBg+Icadeod9dz15K7WNGwgmPzjuXeKfeSYcmI+ryO1lVs2nQVQgQZM/pl0tKOjUG0knTkgrW1uBcvoe399/GuXo2i15N4+mkkX3gh5jFj5AtCSfqR9lbs4os5z9BQvpPCUWM56TdXyQ2VUo8jk+84cLvdXHvttbS2tqLT6RgwYAAvvvgiOp2On//85zz++OPYbDbmzJkT71BjSgjB+xXv8+cVf0YVKvdMvodzB54bk0Sktm4+27ffg9mcz+hRf8di6R+DiCXp8KgeD56VK/EsWYpn8WIClZUA6Pv1I/Om2SSdcw661NT4BilJPZjP42bxvH+y4fOPsCQmcfo1sxky7Tj5QlbqkWTyHQfjxo1j6dKlB73tww8/PMrRHB12n537lt3Hl1VfUppZyp+m/Yl+Cf2iPm847GXnzgeprZtLauqxjBj+JHp9dK0JJemHCFXFt3VrZ7LtXbcOgkEUkwnLxAmkXHQh1mnTMBQXy+RAkqIghGDroq9Y+K+XaW9rY+wpZzLlgp9jstriHZok/Wgy+Za63IKqBdy77F5cARezx83mkmGXRN1CEKDVuYaysptpb6+isOC3FBfPRqORf9JS1wg2NkaS7SVL8CxdSthuB8A4ZAipv7wE27RpmEtL0XS0EJUkKTotNdV88dIz1JRtJnvAIM659Y9kFQ+Id1iSFDWZqUhdJqyG+dvav/HKllcYkjqEOSfPYWBK9K0TVTVAxe4n2bPnBUymXErHvk5KysQYRCxJ31F9Prxr1uBZvATPkiX4d+wAIj23rdOmYps6FeuUKegyot+vIEnSd0LBICvffYMV77yJwWRi5m+vYeTxJ6N0zMCQpJ5OJt9Sl3AH3Pzhmz+wqHYRswbP4pYJt6DXRj/a1+XeRlnZTbjdW8nNuYCBA29Hp0uIQcSSBCGHg7b338e98Bu8q1cj/H4UvR7zuHFk3jQb69SpGAcPlkmAJHWRmq2b+fzFp7HX1TB02nEc98vLsSTJSa1S79Ink28hRK+twxT7jJ2Ol6q2Kq5dcC1VbVXcOelOZg2ZFfU5hQhTVTWH8oq/odcnMmrUi2SknxiDaKW+TgiBb+NGHK/Ppe3jjxGBAIbiYpJnXYBt6lQsEyagkQM7JKlL+TxuFv37FTZ++QmJGVmcc9u99B8zLt5hSVKXiFvyrShKP+A1IBtQgReFEE8oipIKzAeKgErgAiGEQ4lky08ApwNe4NdCiLVH+nNNJhMtLS2kpaX1ugRcCEFLSwsmkyluMayoX8GNX9+Ioii8MPMFJuZEXw7i9e6hbOvNOJ1ryMg4hSGD78dgSItBtFJfpnq9OD/8MDL0pmwrGouF5PPOJXnWhZgG9842n5LUHe1cuZQvX34eb2sr4848m6nn/xx9HJ/HJKmrxXPlOwTMFkKsVRQlAVijKMrnwK+BL4UQDymKcitwK3ALcBowsOMyCXiu4/qI5OfnU1NTQ1NTU4x+je7FZDKRnx+f4UHzts3joZUPUZRYxFMnPEW/xOi6mQghqKubx85dD6IoWoYN+yvZWT/tdS+apKMnUFMb2TC5eDGepUsjQ28GDSL7j/eQeOZP5Bh3STqKfG43C155ga2LviKjqJiz/3C33FAp9QlxS76FEPVAfcfHLkVRtgJ5wE+B4zoOexX4mkjy/VPgNRGpq1iuKEqyoig5Hec5bHq9nv79ZQ/oWHL6nTy44kE+2v0RM/Jn8NCxD2EzRNcGyu/fy9Ztt9HSspDUlKkMHfoQJlNujCKW+oqw24N35cpIsr1kCYE9ewDQ5eSQePppJJ19NuaxY+ULOkk6yirXr+HT55/A42xl8nkXMensWXJCpdRndIu/dEVRioCxwAog69uEWghRryhKZsdheUD1Pt9W0/G1I0q+pdhaUruEu5fcjd1n5+oxV3PFyCuibiO4d++HbNt+N6rqY9Cge8jP+wWKIje4ST9MqCq+LWWdq9ve9eshFEIxmyP9t39+caT/dv/+MuGWpDgI+NpZ+M+X2PjFJ6TlF/Azudot9UFxT74VRbEBbwO/F0K0/Y8nxIPdcMDuQkVRfgv8FqCgoCBWYUrf4w16eWzNY8zfPp+SpBKeOvEphqUNi+qcwWAr27ffw97GD0hMHM2woY9itRbHKGKpN/OsWEnr/PmR/tutrQAYhw0l7dJfY50yBfO4cWgMhjhHKUl9W03ZZj557nGcTY2M/8k5TL3gF+jk/VLqg+KafCuKoieSeP9bCPGfji/v/bacRFGUHKCx4+s1wL5FxPlA3ffPKYR4EXgRYPz48fFv/dELrW9czx2L76DaVc2vhv2Ka0uvxaiNbrBIS8tCyrbeSjBop7j4RgoLficH5kg/yL9zJ42P/hX3woVoU1OxzZiBddpUrJMno0tPj3d4kiQBgXYvi+a+yvpPPyQpK5tZ9/yZ/KEj4h2WJMVNPLudKMBLwFYhxGP73PRf4FfAQx3X7+3z9WsURZlHZKOl80jrvaXoBMNBnt3wLC9vfpkcaw4vn/Iy47PHR3XOUMjDrvKHqK19Hat1IKNH/53EBPmgLP1vwcZGmp96ita3/4PGaiXzptmkXHKJnC4pSd1M5Ya1fPbiU7hamhl72k+YduEvMZjM8Q5LkuIqnkuLU4FLgE2Koqzv+NrtRJLuNxRF+Q1QBZzfcdtHRNoM7iLSavDSoxtu37bDsYPbF93Odsd2zhl4DjePvznqTZWR8fA30d5eTUHB5RT3vxFtlCvoUu8Wdnuwv/wyLf/4ByIUIuUXPyf9yivRpaTEOzRJkvbhc7v5+p9z2PL1F6Tk5nPhH/9C3pDoShMlqbeIZ7eTxRy8jhvggOkpHV1Oru7SoKQDhNUwr5a9ytPrnibRkMjTJzzNjH4zojqnEILqmlfYufNBOR5eOiz+nTtxzJ2H8733UD0eEk47lcwbbsAg93VIUrciVJWdK5ey4B8v4G1zMvFn5zP53Itkbbd01ASDQRRFQdeNu+d038ikuKtuq+aOJXewrnEdMwtnctcxd5Fiim6FUVUDbN9+D3X1b5CefhLDhz0qx8NLByUCAVxffIFj7jy8q1ah6PUknHYqqZdcgnnkyHiHJ0nSPrxtTjZ/9Tkbv/wE594GMgr7c/Yt98hOJlLUwuEwXq8Xj8dzwOVgXw8EApx77rmM7MbPEzL5lg4ghOCtnW/xyKpH0Ck6/nzsnzmj/xlRt2YLBFrYtOlqWp2rKCq8kuLiG2ULQekAwfp6HG+8QeubbxFubkafn0/mTbNJOuccdKmp8Q5PkqQOQgjqd25j/WcfsWP5YsLBIPnDRjDtwl8ycOIU2bdbOoAQAr/fj9frPejl24R63499Pt9Bz6UoClartfOSnJzc+XFmZuZBv6e7kPcMaT9N3ibuXno3i2sXc0zOMdw/9X6yrdlRn9fl3sbGjb8lEGhm+LDHyc4+KwbRSr2FUFU8S5bimDcP91dfgRDYpk/v7MutaOSLNEnqTup3bueLl56lcXc5BrOZkSecwuiZp5HerzDeoUlHmRACn8+H2+3e7+LxeA649ng8hMPhg55Ho9FgsViwWCxYrVZycnKwWq2dn+97sVgsmEwmND30uUEm31KnTyo/4U/L/4Q/5Of2Sbcza/AsNDFYmW5q+owtZbPRaRMYVzqPxMRRMYhW6g1CDgfO/7yDY/58glVVaFNTSbv8cpIvuABDfl68w5Mk6XvCoSDL357HinfexJaaxswrrmHItBmyg0kv5vV6aW5upq2t7aCXQyXUGo0Gq9WKzWbDarWSlZW1X/L8/YvRaOwzw89k8i3h9Dt5YMUDfLz7Y0alj+KBaQ9QlFQU9XnDYR/lFY9SXf0PEhNGMWrU8xiNWdEHLPVoQgh8GzfieH0ubR9/jAgEMI8bR8Z115Fw8kw5DEeSuqnm6j18/MxjNO4uZ/iMEzn+17/FaLHGOyzpMAghCAaD+Hw+fD4fgUDgoMcFg0Gam5tpamrqvLjd7v2O0ev1JCYmkpiYSP/+/bHZbPtdvk24zWZzn0mmj5RMvvu4pXVLuWvxXdh9dq4dey2XjbgMXQyG27S1bWRL2U14veXk5f2CgQNuQ6s1xSBiqadSvV6cH35I69x5+MrK0FgsJJ17DikXXohp8OB4hydJ0iEIVWXNR++xeN5rGExmzpp9OwMnTol3WFIHIQRutxuHw0Fra+t+1263m/b2dnw+3yHLPQ7GYDCQkZHBgAEDyMjIICMjg+TkZBITE/vUCnVXkcl3HyWE4NUtr/LYmscoSS7h6ROfZmja0KjPq6pBKvc8R2Xl0xgMGYwZ8yppqdNiELHUU/krKiJtAt99F9XlwjhwINn33E3iT85Ca5OrZpLUnbXUVPHFS89SU7aZkvHHMPOKq7Emy776sRYMBmlpaaGlpQW73U4oFDrocT6f74AuH16vF1VV9zvOarWSkpJCZmYmJpMJs9m837XBYDhoAq3VaklLSyMxMVEm2F1IJt99kD/s575l9/Hf8v9yStEp3D/1fsy66Ov1PJ5ytpTNxuXaRHbWzxg06G70+qQYRCz1NCIYxPXlAhxz5+JdsQL0ehJPPpmUiy/CXFoqH9QlqZtzO+wse/N1Ni34DL3JxClX/p7hM06U991DCIfDuFyug96mqirt7e0HXDweT2fC3draelg/x2Aw7NfdIy8vD4vFQmJiIsnJyZ0Xgyzf69Zk8t3HNLc3c/1X17OxaSNXj7ma3436XdQPpkKo1NS8xq7yh9FqLYwY8TRZmafFKGKpJwk2NND6xpu0vvkmoaYm9Lm5ZNxwA8nnnoMuPT3e4UmS9AMCvnZWv/8fVr//DuFQkDGnnsEx51yIJVEupHzL5/PR0NCw36WpqemIyjoAjEYjqamp5OfnM2bMGNLS0khPTyctLU0mz72cTL77kLKWMq5bcB1tgTYeO+4xZhbOjPqcPl8dZWU342hdTnraCQwZ8iBGY0YMopV6klBzM03PPEPrG2+CqmI9dhrZ992Lbfp0FK023uFJkvQD1HCYzV99ztI3/42n1cGgSVOZdvGvSMnOjXdoR004HKatrQ2Hw4HD4cDlch10mIvX6+38HqvVSnZ2NiUlJaSmph50MUuj0WAymbBYLJjN5s6LVj429lky+e4jPqn8hLsW30WyKZnXTnuNIalDojqfEIKGhv+wfcd9gGDokD+Tk3O+fEuyj1G9XlpeeQX7nJdQAwFSZl1A6qWXYujXL96hSZJ0mOx1tXzyzGPU79pO7qChnDX7dnIHRb8HqLtSVRW73U5NTQ11dXWdddZOp/OA2mmz2dzZZzo9PZ3CwkKSkpLIzs4mOzsbm80mn/ekIyaT717OHXDzl1V/4d1d7zImYwyPH/846ebo3v4PBJrZtu1Ompo/JzlpAsOGPYLZLJOtvkSEwzjfeYemJ58i1NhIwsyZZNx4A8b+/eMdmiRJh0moKus//4hv/vUPdHo9p197E0Omzug1yeS3w19aW1tpbW2lrq6O2tpaamtr8fv9QKSGOj09ndzcXIYPH05KSgopKSmkpqaSkJAgV6elLiGT715sVcMq7lx8Jw3eBq4YeQVXjr4SvVYf1TlbWr5hS9lNhEIuBgy4jYJ+l6Io8sGpr1D9flyffkrL3+fg37kT85gx5P3tcSylpfEOTZKkI+BqaebT559gz8Z1FI0Zxym/uw5balq8wzpigUCgs0zEbrd3ttlzOp20trbu189aURSysrIYMWIE+fn55OXlkZ6e3mOnJEo9l0y+eyFfyMeT657kn2X/pDCxkNdOe43RGaOjOqcQgurql9m56yFstkGUjv0nNpvszdxXBGpqaJ0/n9a33ibscGDo35+8v/2NhFNO7jWrZJLUFwgh2Lb4a758+XnC4RAnXX41o046tVvfj71eb2di/f3r73cYMRqNJCcnk5KSQlFREcnJySQlJZGcnExGRobcyCh1CzL57mW2NG/h9sW3U+Gs4KIhF/H70t9j0VuiOqeq+tm27S7qG94mI+NUhg97BK02unNK3Z8Ih3F/8w2OefPwfLMIFIWEE08g+cILsU6ejCJXiySpxwiHguxatZz1n35IzdbN5A4ayqlX33BUN1S2t7cfcrLivuUh+14cDgc+n2+/YxMSEkhJSaGkpKSzRCQ1NZWUlBQsFvncJHV/MvnuJYJqkDmb5vDihhdJNafywswXmJIb/QQyf6CZTZuuxOlcS/+i6+jf/1oURSZdvZkQAs+iRTQ+8ij+nTvRZWSQfuWVJF9wPvrs7HiHJ0nSEWhrbmLTl5+w8ctP8TpbScrM4vhf/5Yxp5yBRhP7ksFwOEx9fT1NTU37rVLb7fYDkuhD0el0nf2q8/Ly9kuuU1JS5Oq11OPJ5LsXqHBWcMeiO9jcspkzi8/k1om3kmSMviery7WFDRt/RzDoYMSIp8jKPD0G0Urdma+sjL2PPIJ32XL0BQXkPvooiaecjKKPbq+AJElHV+32raz679tUrFmJQFA8djxjTj6DotGlMX/XyuFwUF5eTnl5ORUVFZ2bGRVF6SwBGTFiBCkpKZhMpoOew2AwkJKSQnJyMlartVuXwUhStGTy3YOpQmXutrk8vuZxzDozf53xV04uOjkm525o+C9bt92OXp/M+HFvkJAwPCbnlbqnYF0dTU88gfO/76NNTCTr9ttIufBCFLnCJEk9SigQYPH8f7Lmw3cxJyQy4afnMurEU0nKzDric3k8HioqKqitrT2gBR9ERqLv2bMHu90OQGJiIsOHD6e4uJicnBySk5NltxBJOgiZfPdQ9e567lpyFysaVjAjfwZ/nPLHqFsIAgSDrWzffg97Gz8gKamUkSOelUNzejF/RQWO1+fS+sYbAKRd/hvSrrgCbWJinCOTJOlI7a3YxcfPPEZLTRWjZ57O9F9cisFkPuzvD4VC1NTUdK5i19XVAaDX69HpDkwXFEUhLy+PiRMnUlJSQnp6ulyxlqTDIJPvHuirqq+4ffHtqELl3in3cvaAs2PygNfSspCyrbcSDNopLr6RwoLfodHIP5HeRgSDuL5cgGPePLzLl4NeT9KZZ5Jx7TXoc/vONDtJ6i3UcJiV777JsrfnYklM4tzb7qVg1FjC4TDBYPCA430+33612N9+3NTURCAQQFEU+vXrx/HHH09JSQm5ubmyHZ/U7YU9QQKVTvy727COz0KfbY13SIckM6seRAjBS5tf4sm1TzIsbRiPzHiEfgnRD7cJhTzsKn+I2trXsVoHMmb0HFlm0gsF9+6l9Y03aX3jDUJNTehzc8m44QaSzzsXXVrP6+8rSRLY62r4+JnHaNi1g+LJx5I7aTqrdlYw/9MvaG9v/8HvVxSFpKQkUlJSGD16NMXFxfTv3/+QtdmS1F2EXQH8u52RS4WT0F5v5AadBkO+TSbfUvR8IR93L72bj3d/zGn9T+O+Kfdh0kX/4NjqXENZ2U20t1dTUHA5xf1vRKs1xiBiqTsQQuBdtgzH3Hm4FiwAVcV67DSy77sX2/TpKLIeU5J6pDZ7C9+8PZ+ylSsQiUkYJsxgQ6uHDZ9+is1mY9CgQYcsA9Hr9Z3dQ5KTkw9aUiJJ3U3I6SdQsU+y3Rx5cakYNBgKE7GMzsBYnIQhPwFF173fqZH3uB5gr2cv1391PVtatnB96fX8ZsRvoi4zUVU/FbufZM+eFzGZcikd+zopKRNjFLEUb2GnE+e77+KYO49AZSXa5GTSLv01ybNmYegX/bslkiTFVjgcprq6murqakKh0AG3CyFoa2ujpaWFxoZ6/MGOY/L6o9VqycnMYuLkEkpKSsjKypK111KPJoQgbPfh393WubodtkdaVSomLcaiJKwTsjEWJ6HPtaJou3ey/X0y+e7mNjVt4vqvrscT9PDE8U9wQsEJUZ/T5d5GWdls3O5t5OZcwMCBd6DT2WIQrRRv7Zs245g3l7YPP0L4fJjHjCH34b+QcMopaIzyHQ1J6i6EENjt9s7Njbt37z7kAJpvmQwGVE8bqsdNVno6o6dOp2DgIDIzM2Xva6lHE0IQamrvTLQDFU7CbZH7g8aiw1CUhG1KLsb+SehzrCianv3iUibf3ZQQgjd3vMlfVv6FDEsGz898nkEpg6I8Z5g9VXOoqHgcvT6J0aP+Tnp69Mm8FF+qz0fbRx/jmDsX36ZNKGYzST/5CSkXX4Rp6NB4hydJfVooFMLpdB6wubGxsZHW1lYAkpOTGTVqFCUlJRQVFWE279+hpKmqkk+fe4K9W3eS2b+EGZdfQcGI0fH4dSQpJtRAmGCNG/+eNgJ72ghUtaF6I+/maBL0GPsnRS7FSegyLD0+2f6+I06+FUVJAfoJITZ2QTwS0Oht5O6ld7OkdgnH5BzDX6b/hVRTalTn9Hr3ULb1JpzOtWRknMqQwfdjMER3Tim+ApWVOObNp/Wdd1CdTgwlJWTdeSdJPz0LbUJCvMOTpD7D7/d3JtXf7yLidDoRQnQeq9frSUlJIScnhylTplBSUkJqaupBy0RUNczq999h6Rv/wmi1cdo1sxk6dUbMh+RIUldTA2ECu534drXi3+0kWOcBNXK/0GWYMQ1Lw1iYiLF/Eto0U68vmzqs5FtRlK+BszqOXw80KYqyUAhxYxfG1id9svsT7l9+P4FwgNsn3c6swbPQRDnOvbl5AZu33ICiKAwf9hhZWWf1+j/s3izY0EDT357A+d57oNWSMPMkUi66CMuECfLfVZK6kMPhoLq6+oBVbI/Hs99xZrOZ1NRU+vXrx6hRozo3N6ampmKz2Q7rftraUM8nzz1O7bYyBk6cwklXXI0lMfrJxZJ0NIiwIFDrwr+zFd+uVgJVbRAWoFUwFCSQMD0fQ2EChoJEtNa+N0FZ2fcV+SEPUpR1QoixiqJcTmTV+x5FUTYKIUZ1fYg/3vjx48Xq1avjHcZhcfqdPLD8AT6u/JhR6aN4YNoDFCUVRXVOIQR7ql6kvPwREhKGM2rkc5hMso9zTxV2uWj5+xzsr74KqkrKL35B2mWXosuQQ5AkqasIIahwOHlv+UpWVdUQ0EbWrEwmExaL5buL2YLFGvlYr//xyYQQgrod29i1chmKRsOgY6aSVTxAvrCWfhQNYNVqsWg1WLQarB3XJo2GWP5Fqf4wwTo3gTo3wXoPwXoPBCNTUXWZZgyFiRgKEtHnWtEatOgUBa2ioFPo/FgDh4xJ2ef2zo+Vjs9ROr7+7efE7f6iKMoaIcT4HzrucMtOdIqi5AAXAHdEFZl0gHWN65j99WwcPgfXjr2Wy0Zchi7K4TbhsJ9t226jYe97ZGWeydChD6HVHv6kM6n7EMEgjvlv0PzMM4QdDhLPPJOM3/8eQ35evEOTpB7FF1ap8gWobPezpz1Alc9P8CDrT6oQNAaC7Pb62e31EUABcyYMzjz4iVXAo4LHDbijD1STCMecAsBnKrCrLvpzStLRoAcKgILv5xtu2OuGvUcnjKeGFnB+dvctrT3cDO9e4FNgsRBilaIoxcDOrgur76hx1XDtgmtJMiTx9BlPMyxtWNTn9Pv3snHTlbS1baC4+EaKh2DxdAAAIABJREFUCq+SqyY9kBAC1+ef0/TXxwjs2YNl0iQyb74Z8wg5AEmSDsUZDFG5T4Jd2e6nsj3AnnY/df4g++bakRXAgz82JoRD6O3NDHE7GZRoY+aoEYzKyiBZH/ve+O2uNlZ/8B82f/UlGq2GyedezMgTZsrabilqISHwhtXOi6fj2qeqR3QeNaQSamonuDeyqh1u9QOg6DXoMszo0syd15of6LGtAmEhCAlBWERiDAmBeohKDAEIEblWiTw3CiIl44LIx/seIxAMs3XvxcbDTb7r9y0xEUJUKIryWBfF1Gd4g16u/+p6VKHy3EnPUZBYEPU529o2snHj/xEKuxg18jkyMk6OQaTS0eZdu47GRx6hfd06DANKyH/+OWwzZsgXUVKfERaCOn+QPR3JcySJjiTUTYED+2AD+FSV1lB4v69lGHQUmYxMTrZRZDZSZDZ0XBtJ02sPuE+1t7fzr3/9i9raWvLz85k5cyaFhYVd8jsGA37Wffw+K999k0B7O+NPmMmU8y7GlionzkrxJYQg1OjFt92Bb4cD/25nZ822sTAR46h0jAOSMeQloGjl89KROtzk+ymg9DC+Jh0mIQR3LbmLXa27eO7E6BNvIVSqa16lvPwRDIZ0xo15kwTbkBhFKx0tgcpKGv/6GK7PP0eXkUH2/feRfPbZKHICndRL1foClLnb2eMLsNvbsUrt81PVHiCwb5cQRaHAZKDQbGBEgpmDra3pFIWCfRLsQpMBm+7wV6qDwSBz586lvr6ec845h5EjR3bJC16hqmxd/DWL5/0TV0sTxaUTOPbiX5Per2uSfEn6IUIIws4AgWoX/h0OfDvshJ2RPtu6LAu2KbmYBqZgKEpEY5CTkaP1P5/RFUWZDEwBMhRF2bezSSIg/+9HYc6mOXy25zNmj5vNlLwpUZ2rvb2WrVv/gKN1OelpJzB06EMYDHLlpCcJ2e00P/MsjvnzUQwG0q+9hrRLL0VjscQ7NEmKKU8ozNJWN1/bXSx0uNjl9XfeZtVqKDIbGGw1cXJaEkVmA/3NRgrNBvJMBrRd+M5POBzmrbfeoqqqivPOO48RI0Z0yc/Zs3E9C//9Mk2VFWQVD+DUq26gYES37l0g9TJCFQT3egl2bo6MXH/bZ1sxajENTMZ0YirGQSnokuWAtlj7oeU0A2DrOG7fxsFtwHldFVRv903NNzy17inOKD6DXw3/1Y8+jxCC+oa32bHjfkAwdMifyck5X5Ym9CBqezv2V1+j5e9/R/X5SD7/PDKuuQZdenq8Q5OkqAghcITCnaUi5V4/S1vdrHJ6CAqBWaNwTLKNS3LTKE20UmQ2kK7XxeXxSwjBBx98wPbt2znttNO6JPFuqqrkm3//g8r1a0jMyOT0a29iyJTpsq5b6nJCCELN7fh3Rdr++cudCF9H6ZZOgz7bgnlEOvocK/pcG4Z8W48b197T/M/kWwixEFioKMorQog9RymmXq3CWcEt39zCkNQh/HHyH3/0E00g0MzWbXfQ3PwFyckTGTb0YczmfjGOVuoqIhzG+d5/aXrySUINDdhOOIHM2TdiLCmJd2iSdNhUIaj3B/ff3Oj7rj67LbT/pq7hNhNX5GdwXGoCE5OsmLrJE/yCBQtYt24d06dPZ9KkSTE9t8vezNI3/s2Wr7/EYDEz4xeXMebUn6CLoh2hJP2QUKsff0Uk0fbvaiXsjLzDpE02YhmZjrE4CX2eDV2aWdZsx8HhFpIaFUV5ESja93uEEHI2+RFwBVxcv+B6DFoDTxz/BCad6Uedp7HpU7Ztu5NQyM2AAbdR0O8ylCgH8UhHj3vxEhoffRT/tm2YRo0i75GHsUyYEO+wJOmgAqpKtS9AZXuA3e3+/TZAVvkC+NXv6rJ1CvQzReqtx3WsZhd1lI0UmIxYukmyva/ly5ezaNEiSktLOf7442N23m+nUy57ay5CDVN6xk+ZdPYFmG1y+qwUeyGHD3+FM3LZ7SRs9wGgsegwliRjHNAP04BktKm9f3pkT3C4yfebwPPAHCD8A8dKByGE4I7Fd1DjqmHOKXPIseUc8TlCIRc7dtxHfcN/sNmGUTr2UWy2wV0QrRRraiCA69NPcbw+l/Z169Dn55P32F9JOO00+UAodSshVbC2zcPXDhdf211scHkJ79MBzKyJ1GUPsJg4KS2xs3NIkdlAntGA7hCt+7qjTZs28cknnzBkyBDOOOOMmN0XWxvq+fiZx6jbsZUBEyZz3C9/Q1JmdkzOLUkQGdfuL2/Ft8OBb7tjv2TbUJSEbUpuZHU724rSg+6TfcXhJt8hIcRzXRpJLzd321y+qv6KP0z4A+Oyxh3x99vtSynb+gcCgUaKiq6mf9E1aDSGLohUiqVATQ2t8+fT+tbbhB0O9IUFZN1xB8mzLkBjkP9+Uvy1BkNUtgfY4PLytd3FYocLV1hFA4xNtHB1v0xKLKbOzY8ZhvjUZceC3+9nz549lJeXs2vXLlpaWigsLOTcc89Fq42+h4AQgo1ffMLCf76ERqvl9GtmM2TacT32/5fUfYhgmGBTe2fC7a+ItP5T9BqMJcmRZLskGX2WRSbbPcDhJt/vK4pyFfAO0Lk1XQhh75Koepnt9u38dfVfmZ4/nV8M/cURfW847KO8/BGqa17BbC5iXOkbJCWN6aJIpVjxLFtGyyuv4PlmESgKthOOJ+Wii7BOniw3WElHlSoEewPBzlKRb2uzd3d8vG9f7Dyjnp9mpjAjNYFjU2wk63t+i8u2tjY2btzIrl27qKqqQlVVdDodRUVFjB8/ntLS0qjGwX/LbW/h0xeepHL9GgpGjuGU/7uexPSMGPwGUl8igmpkRPteD6HGdkJNXoJN7YQdPr6dEKXLtGCbnItpcArGoiQUvXxO6WkUcYiJQvsdpCi7D/JlIYQojn1IsTN+/HixevXquMbgDXq58MMLcQfcvHXWW6SaDn/caVvbRraU3YTXW05+/iUMKLlFjojv5nzbt9P4yKN4Fi9Gm5FOyvnnk3z++ehzjrzMSJIOV1AV1PgiNdmdCbbvu8mOvn3qsrUK5BsN3xs4Y2Cw1Ux/s6FXrNIKIdi9ezerVq1i27ZtCCHIysqipKSEkpISCgoKYpJwQ6S2u+ybr1j42hxCwSDTf3EpY2aeLl9kS4cl7AkS2NOGf08bgco2AjUuOuu8dBr0GWZ0mZbIdYYFQ0ECupQft19M6nqKoqwRQoz/oeMOa1lDCNE/+pD6podXPUyls5IXT37xsBNvVQ1SWfkMlXuexWDIYMyYV0lLndbFkUrRCDY00PTEkzjffRdNYiKZt9xCys8vlqUlUsx4wuH9xqV/m2TvbvdT6w98ry47Mmymv9nAcakJFHV8XGQ2kmc0oO+lb0u3t7ezfv16Vq9eTUtLC2azmSlTpjBu3DhSUw9/4eNwVW5cxzf/epmmPbvJGTiYU6+6kdTcvJj/HKn3UP2hzg4kvl2thBq9kRu0CoY8G7apuRgLE9Hn2NAmG2UJSS91WMm3oigW4EagQAjxW0VRBgKDhRAfdGl0PdwnlZ/w9s63uXzk5RyTc8xhfY/Hs4stZbNxuTaTnfUzBg26B70+sYsjlX6ssNtNy9/nYH/1VQiHSb30UtJ/91u0SUnxDk3qgezBEJVef8cKdmT1+tsE+/sj1VN0WgrNRsYlWjjXnELhPmPTs3pwXfaREELQ0tJCeXk55eXlVFRUEAqFyM/P5+yzz2bYsGExW+HeV9Oe3ZGe3RvWkpiRxRnX3czgycfK1W7pACIUKSPx72zFt9NBoMoFaqRW29A/CcvYTIxFiZHe2no5u7CvONyCvn8Aa4hMuwSoIdIB5agn34qinAo8QWTC5hwhxENHO4bDUeuu5b6l9zEqYxRXjbnqB4//bjz8w2i1VkaOeIbMzFOPQqTSjyGCQRxvvEHzM88StttJPPNMMn7/ewz5ctVLOnzesMqyVjcL7ZHOIju8vv1uzzHqKTIbIl1FTJGWff0tRopMBpJ6QT32jxEIBNi1a1fnpkmn0wlASkoKY8eOpbS0lJwuKvNqa26K9Oz+5ktMFiszLvkNY045U/bsloDI5MhQo5dAjYtAjZtAjYtgvSdSRqKAPtdGwvR8jAOTMRYmoujki7W+6nAfvUuEELMURbkIQAjRrsRhWUVRFC3wDDCTyAuAVYqi/FcIUXa0Y/lfgmqQW765BYHgL8f+Bb3m0A/MQoRpa9vErvKHaW1dQXraCQwZ8iBGo9yo0x0JIXB9/jlNf32MwJ49WCZOJPPmmzGP7JpR1FLvEVIFtf7IKvZmVzsLHS5WtHoICIFJo3BMko0LslMYZDVRaDZSYDJg7oZ9seOlqamJVatWsWHDBvx+PwaDgeLiYqZNm0ZJSUmXlJVA5D5ftWkDGz7/iF2rl6PRaBh/5tlM+tkFmGy2LvmZUs8RavXj29aCb6sd/24nIhAZLKUYtR1lJHkY8m0YS5LRWuWLNCnicJPvgKIoZjr22iqKUsI+XU+OoonALiFERUcc84CfAt0q+X5u/XNsaNrAI9MfIT8h/4Dbfb467PbFtNgXYbcvIRRyotVa5Xj4bs67dh2NjzxC+7p1GAaUkP/8c9hmzJD/XlInb1hlz/enPXr9VPr81PgChPapyx5iNXFZfjrHpSYwKckmE+2DCIVCbNu2jdWrV1NZWYlWq2XYsGGUlpZSUFAQk/aAh+Jzu9my8Es2fP4RjvpaTAmJjD/zbMacfAaJGZld9nOl7k2ogmCtm/atkYQ7WO8BQJtqwjIuC0O/BAz5CejSzbJeWzqkw+12MhO4ExgGfAZMBX4thPi6S6M7MI7zgFOFEJd3fH4JMEkIcc3Bjo9Ht5NAOMDPP3wFl/5gXUkEGm0QjS4Y+UzVEgpYCAfMhANmhCqffKUjoyJQRRiB+sMH9yLfPWrF/8lNoOAzGHBabHhN+9/vjcEAiV43SR43SV4Pie1ukrxuUtwuLIGDr1+oGhDySTtCpwGNBsIq+IMQCO77jx81BQWjToPFoMVi0GExaDEbtIhwiOotmwgF/OQMGsKYk89g0KSp6OQG6j5JCEGwzoN3QxPtGxoJOwOggKEwEfPQNExDU9FlmOVCjBTzbiefK4qyFjiGyLPd9UKI5ihj/DEO9pe930Oxoii/BX4LUFBQcDRi2o9Ba8BhsOHQHWKTpNAQDuhQQ3qE0AJKpHpddhCUohR5IS06/uu9DvX0Fs/f2RgK0q+lgSSvuzPBTvR6MIUCh/we9SCLtsq3FxH5NxR9/bncFwCXC3y+Hz72RxBC4FIFzn1aMaKATqPBlTyYqswxuKxZsApYtaxLYuhKNpOOglQLhakW+qVaKEi1UJBmIcNmlIniYVBb2mnf2Ix3fSOhpnbQKJgGpZB4SjqmwamyjET60f7nyreiKEOEENsURSk92O1CiLVdFtnB45kM/FEIcUrH57d1xPHngx3fHfp8S1JXCrcF8O104NvpwL+zFdUTeVfFWJyEdXIO5mFpKL2xnCHog9YqsFeAYzfYd8O2D6GtBgadCifdC5lD4h3lEQuHVSo3NrN5YS012xxotAolYzMYMSOfnAFJMmHqAkII7J4Ae+xeqlq87GnxssfuwRcM//A3d2NCgLM9yJ4WL/XOdtTe/Io8hgagYRo6jkXPYCKvkHcYYXuKjr05ZmzJZmwmXTd4z61rKAoYtBqMei0mvQajTotRF7nu8jfkFNAqCjqtglajQadR0CgKWo2CRonEpigKCqBRlMjnRK73jV9RFFIseiyGo78p/XBXvn8o+X6xo7XgVwe5WQghTogmyCOlKIoO2AGcCNQSWY+4WAix5WDHy+Rb6kuEKgjWe/Bts+NZ3UDY4UeTYMA6MRvrxGx0ScZ4h9i1gu2w4nlY9BgE3FD6SzjuNkjIjndkP4qjwcOWb+rYuqyeQHuI1FwrI6bnMXhSNgZz3+x0Iv04gZBKXWt75AWG3YvDc+h3ZPoaRRVkOALkNPnJafJh8akIwJ6kZ3eqnrUWhd2BIM3uAM1uP3ZPgLB8JdPtPXbBaM4pPXDPXVeLSfLdHSmKcjrwNyLFGi8LIR441LEy+Zb6KqEKfDsceJbV4dvhAAVMQ9NIODYPY1Ev70HuaYFvHoZVc0BrhCnXRi7GntmZIugPs3P1XjYvrKWpyoXeqGXQpGxGTM8jPb9n/k6SFG/htgDuZXV4VtSjekMoeg3GAcmYh6VhGpKKNuHg9f2qKvCFeva7Iv+LKiIv1vyhMP6giq/j2h9S6ep8URWgCkFYjVxCqiCsqoRUgRCR0kIhvv1YoKrflRuKjlK9b78wsX8qRenWLo33YGKafCuKcjXwbyFEa8fnKcBFQohno460C8nkW5IgZPfhWVGPZ3UDqieEaVgaSacVoc+wxDu0rtVSDl/eC2XvgS0rsgo+9hLQ9sxVYyEEjZUuNn9Tw87VjYSDKjklSYyYkUfJ2Ey0+l5YXiRJMRaoc+NeXIt3QxOoAtPQNKzjszAOSEZjkENupOjEOvleL4QY872vrRNCjI0ixi4nk29J+o4aCONeVItrYQ0iFMY6MYfEkwrQ2np5B4fqVfDZnVC9HNIHw8x7I3XhPbh+2ucJsm1ZPZsX1uJsasecoGfolFyGH5tLYrrcvS1J+xKqwLfNjntxLf4KJ4pBg3V8NrYpuejk/UWKoVgn3xuB0aLj4I5hNxuFEMOjjrQLyeRbkg4UdgVo+7IKz8p6FJ2WhBn52I7N692rPkLAtg/giz9Cyy4onAYn3w95B91L3mMIVVCzzcGmhTVUbmxGAIUj0hgxPY+C4WloZMtCqQ8LtwXwrG7As7KBcKsfbZIB25Q8rBOy0FhkpxIp9mKdfD8CFAHPE6mo+T+gWggxO8o4u5RMviXp0IJNXpwfV+Ira0Fj0WEZn41tUja6tF68EhQOwppX4OuHwNsME38HJ/0RDD2/BMdl91G2uI6yxXV42wIkpJkYMT2PoVNyMB+iflWSehuhCvzlrXhW1NNeZgdVYByQjHViNubhvbT7k9RtxDr51gC/I9JlRCEyaGeOEKJb7zqQybck/TB/pRP34lray1pAgGlQCtZJOZiGpPbeCW2+NvjqQVjxHKQNgLNfhPxx8Y4qJsJhlYp1TWz5ppbaHa1odAoDSjMZMSOf7OJE2a5Q6pVCrT68axvxrtlLqMXXsaCQhXViDnpZWiIdJb2228mRkMm3JB2+sNOPe2XkLVrVFUCbbMQ6KQfrhKzeWxdesRDevQpc9XDsbJh+M+h6z+9qr/OweVEt25fVE/CFScu3MWJ6HoMmZmEw9cyNp5L0LdUfpn1zM941e/HvdoIAQ/9EbBNzMI9IR5GbkKWjLFZ9vt8QQlygKMomDjJATggxKrowu5ZMviXpyImwSnuZHc+yOvwVTtAqmEekYzsmB0NRL1w59Tnh41thw+uQPQrOeREyh8Y7qpgK+ELsXLWXTQtraalxozdpGTIpm+Ez8kjLle0KpZ5DCEFgdxueVQ20b25GBFW0aSasYzOxlGahSzXFO0SpD4tV8p0rhKhTFKXwYLcLIfZEEWOXk8m3JEUn2OiNtClcsxfhC6PLsmCbnINlbCYaYy9bOd36Prz/e/C7YOr1MPU6MCbEO6qYEkKwd3cbmxbWsGtNI2pIkDswmREz8igek4FWJ1cKpe4p7AniXduIZ2U9oaZ2FKMWy+gMLKWZGAp74aKA1CPFKvleK4QoVRTln0KIS2Ia4VEgk29Jig01EKZ9QxPuZXUE6zwoRi0Jx+Vjm9rLuqS4m+CTW2Dz22DNhONuhdJf9dje4P9LuyvA1qX1bFlUS1uzD6NFR/6QVAqGpdJvWCoJcgVRijMhBIHKNjwr6vFuboaQwFCQgHVSDuaR6b3rsUfqFWKVfG8GHgHuBm7+/u1CiP9EE2RXk8m3JMWWEIJAtQvXV9X4ttrRJhpIPLkIS2lm79qcWbMaPrsLqpZC2sBIb/DBp/fo3uCHIlRBVZmdXWsbqS6z42n1A5CSbaHf0FSKRqaTPzRFrixKR43qC+Fd24h7eT2hRi+KSYtlbCa2STnos4/+1EJJOlyxSr6nAT8HLgD++72bhRDisqii7GIy+ZakruOvaKX1o90Ea9zos60knd4f06CUeIcVO0LA9o/g83ugZScUTIGT/9RruqIcjBACe72H6jI71Vvt1O1oJRRUyS5OYup5A8guTop3iFIvFqhz41lej3d9IyKgou+XgG1SNuZRGXKVW+oRYpV8ny+EeFNRlN8KIV6MaYRHgUy+JalrCVXQvqkZ56eVhO0+jAOSSZiRj7EkufeshIeDsPY1+PrP4GmC4efAiXdDav94R9blQsEw25c3sPL93XjbApSMzeCYn5WQnNXz+6JL3UPI7sO3w4F37V4CVS4UvQbz6IzIBu/83rXnQur9Yl3zvVYI0eNGwcnkW5KODhFScS+vx/VVNaoniC7NhPWYHKzjetEkOb8LljwJy56OJOQTr4i0JrSkxjuyLhfwhVj/RTXrPq9CDaoMn5HHhNOL5PAe6YipvhD+8lZ8O1vx73QQavEBoEs3Rx4zSjN7z2OG1OfEKvn+HNABY4BF379dCHFWNEF2NZl8S9LRJUIq7ZuacS+vJ7CnDXQaLN+uYvXrJatYbfXw9YOw7l9gSIDpsyOTMvW9f4Oix+ln1Qe7KVtSj06vod+wVPoNjWzSTJSDTKRDECGV9s3NeFY1RPpxq6AYNBiLkzENTMY4KAVdulnuK5B6vFgl3wagFPgncPn3bxdCLIwmyK4mk29Jip9AnTvSpWBdR/1mvg3bpBzMo3tJ/ebeMvjiHtj5GST1gxPugpHng6b3t+uz13vY8EUVVWV23I7IBs2kTDMFQ1PpNzyNgqGpaOWAkz4v2OTFs6IB79q9qN4Q2lQTllEZmAYlYyhIRJGtLaVeJtbj5TOEEE2KoliFEJ6YRHgUyORbkuJP9YXwruvoXLDXi2LSYR2fhXVSNvqMXlA7XLEQPr8L6jdAzmiYeT8Uz4h3VEeFEILWvV6qyuxUl9mp3eEgFFAxJ+gZOjWX4dNy5Yp4H/PtKrd7RQOB3U7QKJiHp2GdmN279oJI0kHEOvmeDLwE2IQQBYqijAZ+J4S4KvpQu45MviWp+/h2Mp17eR3tm1tAFZ0bNE0De3iXFFWFzW/Bl/eBsxoGzISZ90HWsHhHdlSFgyo12x1sWVRL5cZmBFA0Io3h0/MoGJ6GRiZevVbI7osM5Fq9F9UTRJtqwjoxG+u4LLRyb4DUR8Q6+V4BnAf8VwgxtuNrm4UQI6KOtAvJ5FuSuqewK4BnVQOeFfWEnQGMg1JIOq0/hpwe3sM36IOVL8A3f4WAC8b8HI6/AxJz4h3ZUeey+yhbXMeWxXW0twVITDcx/Ng8hk7JkRs1ewmhCnw7HHiW1+PbbgfANDQN2+Qcucot9UkxT76FEJMURVm3T/K9QQgxOgaxdhmZfEtS9yZCKu6ldbQtqEb4Q1hKs0g6uRBtkjHeoUXHa4dvHoGVfwetHiZf0yvH1R+OcEilYn0TmxfWUrezFY1OYUBpJiNm5JNdLMeC9zSRbiVOfDsd+LbZ/7+9+46Ps7rzPf45M+q92pJVLMvIvVuyDabaxtiGQCBAyGYDgSSkk2RvlsCyd29ustmEzSbZbJLN3RQSEpLQIRCq6Rjc5CLbuMpdltxkFatrZs7943kMiiPJtizpmZG+79drmNGZZ2Z+OhrGXx2d5xyC9e34kqNJLMshcU4uUWkR/v+uyHno7/D9OPBD4KfAPOAuoNRae8v5FjqQFL5FIkOopZPG1w/S9G41xmdIujiP5Mvy8cVF+LbuJ/Y6U1HeexISs7tsVz88l1I7Ud3MlrcOsWNVDR1tQTLzkphyWR7j5owkJtJ/1kOUtZbOqibadtbRtquOjgMnIWTfX60kYdYI4idnYvw6eVKkv8N3FvBjYBHgA14CvmKtrT3fQgeSwrdIZAmcaKPhpX20VhzDxPhJmJlN0oWjIn9L6ap18PI/D4vt6s9GR1uAXWuPsPnNQ9RWNREd52fC3BwmX5ZH5qgkr8sT3BOl1x2haXUNgaOtAETnJRFXkq7VSkR60K/hO1IpfItEpo5DTTS9W01LxTEIhIgpSiFpXi7xU7Ii9x98a2HHC87yhMd3utvVfxvyz/g5PWRZazmyt5Etbx5i17ojhAKWUSVpTLksj+IZ2fgj9WcdwToOdVkitNPd4n1ODnETM/Anaa6+SG/6e+Q7H/gJMB+wwAqcke+q8y10ICl8i0S2YHMnLeuO0Ly6hkBtG77EaBLLRjpzSzMidFObYAA2/A5e/y40H4XJ17vb1Rd7XZmnWk92sO3dGt57+xCNx9uIT4lh0vxcJl+SR3Kk/qwjhO0M0bLpGM2ra7TFu8h56O/wvRz4I85mOwB/D3zcWnvleVU5wBS+RYYGG7K0V9bTtKqGtm3ObLe48RkkXphLXEl6ZK6q0H4S3v2Jcwl2Qtmn4bK7h8V29b2xIcuBrSfY8mYV+7bUYoDRU7OYelkeBRMzIvNnHaYCta00ra6hpdzZBEdbvIucn/4O3xuttTPO1BZuFL5Fhp5AfTvNa2poXnuY0ElnPeGkuTkklObgT4zAwPBX29UnOcsTln0Kskq8rsxzjbWtbH27mq3vVNN6spOU7HimuMsVxiVF4M86DNigpW3HCZpW1dC+sw58ED8pk8R5o4gdm6rVZ0TOQ3+H71eA3wJ/cps+BtxurV14PkUONIVvkaHLBkO0vldL86oa2vc0QJQhYWo2ifNyiSlMjrwQcWQrvP0fsPUZCHXCmEud0fDxy4bt6iinBAMh9mw4xuY3q6ipbCAq1s/MKwuZsahAq6SchUBdG2276mjfWUdbZQO2LYAvJYakOTkkluVE/tKeImGiv8N3Ic4ygxfizPl+F7jLWnvgfAsdSArfIsND55FmmlcfpnndEWx7kOjcRBLn5ZIwYwS+WL/X5Z2bpqOw4fdQ/htnt8ykHJh9G0y8FkZOHrYrpJxSe6glVWJBAAAgAElEQVSJtc/tZff6YySkxDDnQ2OYeFEuPi119z5rrbMW99Za2nbVETjmrFbiT40htiSd+IkZxE3I0PKAIv2sv8P3g8BXrbV17tcZwH9Ya+8470oHkMK3yPASag/SsvEozatq6KxpxsT6SZw9kuTLC/CnRNhKDaEg7HoZ1v4aKl8BLCSNhLELnEvx5ZA0wuMivXN4TwPvPF7J4T0NpOcmctH1Yxk9NTPy/uLRj0ItnTSvd97/geOtmGgfscWpxJakEzcunajs+GHdPyIDrb/D9/s7W/bWFm4UvkWGJ2stHQdO0ryqhpZNx5yNey7NJ/nSPHyxEThNobEadr8Ou191rludrbzJmeqG8YVQOA+ihtf0AWstezce592nKmk42sqokjSmLyygaGrmsBoJ76g6SdOqGlorjmE7Q8QUJjt/+ZmahYmOsL/8iESw/g7fFcDlp418v2mtnXrelQ4ghW8RCdS2Ohv3bDqOLymalEWjSSzLwfgjdAQwFIKajbD7NedycDWEAhCdAEUXfxDGs0qGzRSVYDDE1rerWffifprr20lKj2XyJaOYOH8UiUN0PrPtDNJScZymVdV0VjVhon0kzBxB4txcYvK0UZGIF/o7fN8K3As8jjPn+2bgO9ba3/f6QI8pfIvIKe0HGml4fi8d+xqJyo4ndckY4iZlRP6f4dtPwr4VThCvfBVO7Hbasyc6q6ZM+yjEpXhb4yAJBUPs21TLlreqOLitDp/PUDwzmymX5TGqJC3yf9ZA5/FWmlfVOOc3tAaIGhFP0txcEmaPxKeTT0U81e87XBpjJgELAAO8aq3den4lDjyFbxHpylpL29YTNLy4l8CxVqKyT61rPBJf/BAJLnX7nDni63/vjJDHJMG0m6H0U5AzxevqBk39kRa2vHWI7StraG8JkJ6byJRL8xg/L4fYCPtZh9oCtO2qo3nNYdp31YPPED85k8R5ucQWa3lAkXCh7eVR+BaR7tmgpaXCOTHt1I5+CTNGOMsUDqU/2R9aB2sfgC2PQ6ANCubBjL+DCxZCar7X1Q2KQEeQXeVH2PLmIY7uP0lUrJ9xc0Yy5dI8sgvCc/dGG7J0HmqibVcdbTvr6DhwEkIWf0oMiXNzneUBI+0EYpFhQOEbhW8RObOOQ03OiZkbjzonqxUkEz81y1kdYmTC0BhVbDkBG/8I5Q98MC0la7wzP/yChTD6IohJ9LbGQXB0fyNb3jzEzrVHCHaGyClO4eKbxzGyKDym5XTUNHPyzYO076wj1BIAIDovibiSdGJL0ogtSo3ccxVEhgGFbxS+ReTshVoDNK8/Qsvaw3QebgHAlxJDXEk6cSVpxJakR+YOml1ZC0e3uSdrvgr733VGxP0xzqY+s2+HcUvAH1nTMs5VW3MnO1YdZuMrB2g92cmCWycwbk6OZ/UEGtppfHk/LeuPYGKjiJ+UQdy4dGIvSMOfpBFukUih8I3Ct4j0TaC+nfZddc6ugJX1ziikgbgJGaQuHUP0iASvS+wfna1wYKVzouZ7T0HjIUjJg9mfhFm3QrJ3gXQwtDZ18OL/bKF6Vz2zloxm3rXFGN/gjSyH2gKcfKOKkysOgbUkzR9FyuUF+BIi/Jc8kWFK4RuFbxE5f6fm37Zuq6XpnWpsZ5DEshxSFo3GnzyERiWDAdj1Eqz9lTMy7ouCCdc4QXz0fIgaQt9rF8FAiLce2cnWt6spmpbFlXdMGvAt60NtAVrWHaHxtQOEmgPEz8gmdXERURlxA/q6IjKwFL5R+BaR/hVs6uDkawdpWlWDiTIkX5pP0qX5+GKG2EYmtbud+eEbHoK2emfFlDGXfrC7ZkbxkFpD3FrL5jcOseKxXaTnJHD1F6aRkhXff88fsnRUnaR9Zx1tu+rpONgIIYgtTiV12Rhi8sPzxE8ROTcK3yh8i8jA6DzeSuNL+2jdfBxfcjQplxcMzXWWO1udZQtPrSFev99pTxsNFyyCmR+HvNne1tiPDm47wUu/3IIxhsWfmkzBpIw+P5cNWtq21dJScYy2XfXYNmfq0qkTKOPGpxMzOmVonNArIoDCN6DwLSIDq31/Iw0vOBv3mJguOwyOGkLLFZ5iLZzY456s+TrseQM6myF3BpR9GqZ8BGIify58/ZEWnv/5JuoOtzD1inwuvH4s0efwl41gYzvNaw7TvOYwwcYOfMkxxI1Pd1YsuSAt8k/aFZEeKXyj8C0ig6Pj4EmaVtXQUnEMAiFiRqeQOC+XhClZmGif1+UNjLZG2PQIrP01HNsGcakw4+NQeoeztX0E6+wIsuqp3Wx6vYq0kQks+uQkRo7peTlCay3tuxtoXlVN69ZasBBbkk7SvFziJmQM6kmcIuIdhW8UvkVkcIVaOmled5Tm1TUEjrfiS4wioTSHpLm5Q/dkOmudFVPW/gq2PgOhThgx6YP54aMvguj+mz89mA5uP8FrD26juaGD2UtGU3p1EX7/B79MhVoDNK874vy8j7XiSzj1884hKjMyv2cR6TuFbxS+RcQbNmRp311P06oa2rY5I6Fx49JJnJdL3PghPBLadBQqHnbmiR9YCcEOiIpzAvjYBTB2IYyYGFEna7a3Bnj7kZ3sWHWY7MJkFn1yEkmhEE2ramitOOZszFSYTOLcXBKmZQ/dv3SIyBkpfKPwLSLeCzR8MAc4dLIDf1osiWU5xE3IIDo3cegG8Y5mZxOfUydrHt/htCfnfjAqXnw5JGZ5WeVZ27P+KNv+tIPRxpLmM5joLnP884bgHH8ROWdhHb6NMd8HPgR0ALuB26219e599wKfAoLAXdbal9z2JcCPAT/wK2vt9870OgrfIhIubDBE69ZamlfW0L6nAQBfYjSxJWnuLprp+FOG5lraADRUOSdq7n7VOVmztQ4wkDsdJn7I2dQnaYTXVXarrbKehhf20nmoiRafobI5QM6iQmZeM0arlYjI+8I9fC8GXrPWBowx9wNYa79hjJkE/AmYA4wCXgHGuQ/bCVwJVAFrgY9Za7f29joK3yISjoKNHbRV1r2/7nOouROA6JwEYse5K2MUpQ7dKQyhIFRvdEfFl8PB1eCLhknXQumnnGkqYRBqO4800/D8Xtp21OFPiyXlqiJiJmbw+kPb2VV+lJKykSz4xASihto67yLSJ2Edvv+qAGOuB2601n7cHfXGWvtd976XgG+6h37TWnuV2/5Xx/VE4VtEwp0NWTprmp2t7HfW0b6/EYIWonzEFqcSV5JG7Ng0fAndryHuT4yJ/JB+bKezqc/GP0J7A2RPhLJPwbSbnVVUBlmwsYPG5ftpLj+MifWTckUBSRflvd/P1lrWvbif1X/ew4jRySz93DSS0mMHvU4RCS+RFL6fBR6x1j5kjPkpsMpa+5B736+BF9xDl1hrP+22fwKYa639UjfPdydwJ0BhYeHs/fv3D8a3ISLSL0IdQdr3NNC+q462XXUEjrb2/gAD/pRYojLj8GfEEZUZT1TmB9cRtfFPRzNsecJZOaWmwtnivmDuB3PEc2eAb2B/0eg4eJJjD2zBdgRJmpdL8oLCHtfm3rPxGMt/s5WYOD/LPjet1+UIRWTo8zx8G2NeAXK6ues+a+2f3WPuA0qBG6y11hjzM2DlaeH7ecAHXHVa+J5jrf1ybzVo5FtEIl2gvo2OfY3YztDf3mmdTV0CtW0ETrQRqG0l1NT5V4f4EqLwnwrkGXH4EqPpbkKHLyHa2QQmOQzmnVsLh9bD9medkzUPb3La4zNg7BUfhPGUUf36su37Gjj+m/fwJUaTdftkorPPvGlQ7aEmnvvvTbQ0dLDg1gmMm9PdP3siMhycbfgesCERa+2i3u43xtwGXAMstB/8BlAFFHQ5LB+odm/31C4iMmRFpcURNePs1wgPtQcI1LYRPNHmhPLaVgIn2ujY30hrxTE4w3hLdE6iO+88zbt558ZA/mznsuib0HTMOUlz96vOPPEtTzjHZU+ECxY6gXz0/PNaT7ytsp7aB9/DnxpL1memEpV6dtNIMvOSuOmeUl78xRaWP7CV2upm5l1bPHRXsRGR8+bVCZdLgB8Cl1lrj3Vpnwz8kQ9OuHwVKAEMzgmXC4FDOCdc/p219r3eXkcj3yIiH7CBEKH2YLf3Bevb/2beuYn2EZ2XhInqJoAbiEo/NdUljqiMQZrmYi0cec/d5v5V2L8Sgu3gj+2ynvgCGDn5rE/abN1+gtqHthKVGU/2p6f2afQ/GAjx1sM72bqimqJpWVx5xyRiImnKj4icN8+nnfT6osZUArFArdu0ylr7Ofe++4A7gADwVWvtC277MuA/cZYafMBa+50zvY7Ct4jIuXt/3vnOOjqqm7odLbfBEMG69vdXajnFlxj1fhD3Z8YTdSqcZ8bjS4ru/6X5OlrgwLtQ6YbxY9ud9qSRXdYTvwKSsrt9eOuW49T+aTvROYlk3TGlx/ndZ8Nay+Y3qljxWCXpOQlc/YVppGRpp0uR4SKsw/dgUfgWERlYobaAO9+8jeCJ1g/mnx9vJdjQ/lfB3cT4iMqIx5952mh5fD+OEDcdg0PlULUWqsqd1VMAMksgvwwK5sDIKeCPpuNQE/V/riQmP5ms26f0Wx0Ht53gpV9uwRjDkjunkDc+vV+eV0TCm8I3Ct8iIl6ygRCBOjeYu3PPnXDu3Cbg/b8/saMg82Oj8WXm9+tKKvVHWnjuvzfRcKyV2UtGU7qsCH9303dEZMhQ+EbhW0QkXNmQJdjYQaC2FdvDPPR+F2h15ovXbILDmzAnq4j1bcKYgDNnPH00pI+BjGLIGOPeHgNpoyHq3OeBt7cGWPHITravOkx2YTILPzmRzFHail5kqFL4RuFbRER60VAFx3dB3V44sQdO7IW6fc51Z/MHxxkfpORDRhHEp0O3izX2bM+xQl7ffiGdwWjmFa9nesHWcNjA8+z5oiAl1/klJLUA0gohrQBik72uTCSsKHyj8C0iIn1gLTQddUP53r++bmvs01O2BJJ4/dAN7Ds5iVEJe1iY/xgpMXX9XPgACXZAY7WzqkxXcanOXwyk/8UkQGyK08exKRCX4lz7+35C8IDxx0BUHETFdrnEOb+0nSvjA58fjN+59kU5bT39wmvc/xjjXvuc21njIXlk37+nPlL4RuFbRETCh7WW7SsP8/ajOwG4+vPTIudkzFAImo9B/QFoOOBeH4JQwOvKhiDr7Pba1gjtjdB+0r3dAMFw628LwU4IdZ750MF0/f/A9FsG/WUVvlH4FhGR8NN4vJW//GwTjcdbWfa5qRROzvS6JJHzEwpCoB0CbR9c22525e2NtWCDznP91XVvz2Od17HWve1+nTVOI99eUfgWEZFw1Hqyg2f+ayMnapq56tNTKJ7R/TrkIhI5zjZ8a90jERGRQRafHMN1X51JVn4yL/1iC7vKj3hdkogMEoVvERERD8QlRnPdV2YwsjiF5b9+j+0ra7wuSUQGgcK3iIiIR2Lio/jQl2eQNz6dVx/cxpa3DnldkogMMIVvERERD0XH+rn6i9MYPTWTN/+4g7ce2UlnxyBtPCQig07hW0RExGNR0X6WfnYq067IZ/PrVTz6nbUc2du3NcVFJLwpfIuIiIQBf5SPSz46juu+OoNAR5Anvr+O1c/sIRg8xyXbRCSsKXyLiIiEkfwJGdzyL3MZP2ck5c/v44n713GiuvnMDxSRiKDwLSIiEmZi46NY+MlJLP3sVE6eaOPRf1vLe2/rZEyRoSDK6wJERESke8Uzs8kZm8qrv93KG3/YQe2hZi6+6QJ8fo2diUQq/d8rIiISxhJSYrj6S9OZsaiAzW9U8exPKmhr7vS6LBHpI4VvERGRMOfzGebfWMKCWydSXVnP498r50SN5oGLRCKFbxERkQgx8aJcPvy1WXS0BXji/nL2b6n1uiQROUcK3yIiIhEkd2wqN91bRkp2PM/9rIINyw9grfW6LBE5SwrfIiIiESY5I44bvj6b4hnZvPtEJa/9bhvBTq0HLhIJFL5FREQiUHSsn6s+M4Wyq4vYvvIwT/9oPc0N7V6XJSJnoPAtIiISoYzPMOdDxVz1mSkcP9jE498r59iBk16XJSK9UPgWERGJcBfMHsEN/zgbgCe/v47KdUc9rkhEeqLwLSIiMgRkFyZz071lZBUk8dIvt7Dy6d0EA5oHLhJuFL5FRESGiISUGD78tVlMmp/L+hf38/j95RyvavK6LBHpQuFbRERkCPFH+7jiExNZ+rmpNNe389h311L+/D5CQY2Ci4SDKK8LEBERkf5XPCOb3AtSeetPO1n9zB72Vhxj4ScnkZGb6HVpIsOaRr5FRESGqPikGK76zBQWf3oyDcdbefQ7a6l49aA25RHxkEa+RUREhriS0pGMKknjjT/sYMVjuzi6v5Er/n4CUTF+r0sTGXY08i0iIjIMJKbGsuzzU5l7bTE71xzhqR+sp6lOm/KIDDaFbxERkWHCGEPpsiKWfm4qdYdbeOx7azm8t8HrskSGFYVvERGRYaZ4RjYfuXs2UdE+nv7BBravqvG6JJFhQ+FbRERkGMrMS+Kme8rIGZvKq7/dxorHdxEK6URMkYGm8C0iIjJMxSVF86G7pjP1inwqXjnIcz+roL014HVZIkOawreIiMgw5vf7uPSj47j84+Op2lbH498rp/5Ii9dliQxZCt8iIiLC5EvyuO5rM2hr7uTx+8s5uPWE1yWJDEkK3yIiIgLAqJJ0brqnlKT0WJ79aQUVr2lDHpH+pvAtIiIi70vJiueGf5xN0dRMVjy6izce2k6gI+h1WSJDhsK3iIiI/JWYuCiWfnYqs5eMZus7NTzynbXUVNZ7XZbIkKDwLSIiIn/D+AzzPjyWa++aQbAzxJM/WM/bj+yks12j4CLnQ+FbREREelQwKYNb/mUOUy/LZ9PrVTz87dVUbdfJmCJ95Wn4NsZ83RhjjTFZ7tfGGPNfxphKY8wmY8ysLsfeZozZ5V5u865qERGR4SUmLopLbxnH9f9rJsYY/vyfG3njD9vpaNOa4CLnyrPwbYwpAK4EDnRpXgqUuJc7gZ+7x2YA/weYC8wB/o8xJn1QCxYRERnmRpWk89H/PYcZVxaydUU1T/z7OhqOaU1wkXPh5cj3j4C7ga5rGF0H/M46VgFpxphc4CpgubX2hLW2DlgOLBn0ikVERIa56Bg/8z9yAR/68gya69t57LtaE1zkXHgSvo0x1wKHrLUVp92VBxzs8nWV29ZTu4iIiHigYFIGN91b5qwJ/pONbHj5gNYEFzkLUQP1xMaYV4Ccbu66D/gnYHF3D+umzfbS3t3r3okzZYXCwsKzqlVERETOXWp2PB+5u5RXH9zGu09WcuzgSa74xASiY/xelyYStgYsfFtrF3XXboyZCowBKowxAPnAemPMHJwR7YIuh+cD1W775ae1v9HD6/4C+AVAaWmpfgUXEREZQNGxfq76zGTWv5TEqj/voe5wM0s/N5WUzHivSxMJS4M+7cRau9laO8JaW2StLcIJ1rOstYeBZ4Bb3VVP5gEN1toa4CVgsTEm3T3RcrHbJiIiIh4zxjB7SRHXfHE6jcfbePx75VTv0qY8It0Jt3W+nwf2AJXAL4EvAFhrTwDfBta6l2+5bSIiIhImRk/J5KZ7SolNiObP/7mBrSuqvS5JJOyYoXxyRGlpqS0vL/e6DBERkWGlvaWTl3/9HgfeO8HUK/K5+MYL8PnDbbxPpH8ZY9ZZa0vPdJz+TxAREZF+FZsQzdVfnM70RQVsfr2KZ39SQVtzp9dliYQFhW8RERHpdz6f4eIbS1hw60SqK+t57HvlnKhu9rosEc8pfIuIiMiAmXhRLtf/wyw624M8dn85u9Ye8bokEU8pfIuIiMiAyilO5eZ7y8jOT+LlX7/HW3/aQbAz5HVZIp5Q+BYREZEBl5Qey3X/MJMZiwrY/OYhnvzBehprW70uS2TQKXyLiIjIoPD7fcy/sYQln51C/eFmHv23tezfUut1WSKDSuFbREREBtXYmSO46d4yktLi+MvPKlj19G6CAU1DkeFB4VtEREQGXdrIBG78xmwmXpjLuhf38/j95dRWN3ldlsiAU/gWERERT0TF+Flw60SWfnYqTXXtPPZv5WxYfoBQaOhuACgS5XUBIiIiMrwVz8wmZ2wqb/xhO+8+Ucm+TcdZeNtEUrLivS5NpN9p5FtEREQ8l5ASw9LPTWXBrRM5dvAkD397DdverfG6LJF+p5FvERERCQvGGCZelEve+DRee3Abr/1uGydPtFF2dRHGGK/LE+kXGvkWERGRsJKSGc+1X5nBhAtzWPuXvax4dBdW88BliNDIt4iIiIQdn9/Hgk9MJDYhmopXD9LW0smCWyfi92vcUCKbwreIiIiEJeMzzL/xAuISo1n9zB46WoNc9enJRMX4vS5NpM/066OIiIiELWMMpcuKuOxj49i3+TjP/qSCjtaA12WJ9JnCt4iIiIS9KZflc+Udkzi8u4Gnfrieprp2r0sS6ROFbxEREYkI48pyWPaFaTQcbeWx767l8J4Gr0sSOWcK3yIiIhIxRk/J5CPfmE1UrJ+nfriere9Ue12SyDlR+BYREZGIkjkqiZvuKSWvJI3Xf7+dtx7eSTAY8roskbOi8C0iIiIRJy4xmmu+NJ0ZiwrY/EYVz/54I60nO7wuS+SMFL5FREQkIvn8PubfWMKiT07k8J5GHvtuOccOnPS6LJFeKXyLiIhIRBs/L5frvz4Lay1P/Ps6tr2reeASvhS+RUREJOKNLErh5n8qI2dsKq/9bjuvP7SdQGfQ67JE/obCt4iIiAwJ8ckxXHvXdGZdNZqtK6p56j/W01jb6nVZIn9F4VtERESGDJ/fx4XXj2Xp56ZSf6SFR/9tLQe21npdlsj7FL5FRERkyCmekc1N95aRmBrLsz+pYO1ze7Eh63VZIgrfIiIiMjSljUzgxm+UMq5sJGue3cuzP63QcoTiOYVvERERGbKiY/0sun0Sl398PNU763nkO2uprqz3uiwZxhS+RUREZEgzxjD5kjw+cvds/NE+nv7hBja8fABrNQ1FBp/Ct4iIiAwL2YXJ3PxPZRRPz+LdJyt5/uebaWvu9LosGWYUvkVERGTYiI2P4qo7p3DxTSUc2FLL4/eXU3e42euyZBhR+BYREZFhxRjD9IUFfPgfZtLRGuDx+9dpOUIZNArfIiIiMizlXpDGjd8oJTkjlr/8pIKK1w5qHrgMOIVvERERGbZSsuK54R9nUzQtixWP7uKNP+wgGAh5XZYMYQrfIiIiMqzFxEWx9LNTmb3E2Zb+mR9vpLVJ64HLwFD4FhERkWHP+AzzPjyWK++YxJG9jTzy7TXs36J54NL/FL5FREREXOPm5PCRu2cTmxjNX35awWu/30Z7a8DrsmQIUfgWERER6SK7MJmb7y1j1lWj2f5uDQ9/azUHt53wuiwZIhS+RURERE7jj/Zx4fVjueEfZxMV4+eZH2/kzT/uoKNNo+ByfqK8LkBEREQkXOUUp/LR+8pY9cweKl49yO4NR5kwL5dJF48ibWSC1+VJBFL4FhEREelFVIyfi28sYezMEWx4eT8bXz3IhuUHGFWSxqSLRzF2VjZR0X6vy5QIYbxaTN4Y82XgS0AAeM5ae7fbfi/wKSAI3GWtfcltXwL8GPADv7LWfu9Mr1FaWmrLy8sH6DsQERGR4ai5vp1tK2vY9k41jcfbiE2IYuL8UZQuHU1sQrTX5YlHjDHrrLWlZzzOi/BtjLkCuA+42lrbbowZYa09aoyZBPwJmAOMAl4BxrkP2wlcCVQBa4GPWWu39vY6Ct8iIiIyUGzIUrWzjq0rqtm97ihxSdHMv7GEcXNGYozxujwZZGcbvr2advJ54HvW2nYAa+1Rt/064GG3fa8xphIniANUWmv3ABhjHnaP7TV8i4iIiAwU4zMUTMigYEIGxxaf5I0/7uCV32xl64pqLv3YODJHJXldooQhr1Y7GQdcYoxZbYx50xhT5rbnAQe7HFfltvXULiIiIuK57MJkbrx7Npd/fDy1h5p49F/X8u6TlVodRf7GgI18G2NeAXK6ues+93XTgXlAGfCoMaYY6O5vNJbuf0nodr6MMeZO4E6AwsLCcy9cREREpA+MzzD5kjyKZ2Sz8qndbHj5ALvWHmHeh8cyrmwkxqepKDKA4dtau6in+4wxnweetM6E8zXGmBCQhTOiXdDl0Hyg2r3dU/vpr/sL4BfgzPnu8zcgIiIi0gfxyTEsuHUiE+eP4u1HdvLKb7ZS8epB5t94AXnj0r0uTzzm1bSTp4EFAMaYcUAMcBx4BrjFGBNrjBkDlABrcE6wLDHGjDHGxAC3uMeKiIiIhKXcsancdE8pi26fROvJDp7+4Qae++9N1B1u9ro08ZBXJ1w+ADxgjNkCdAC3uaPg7xljHsU5kTIAfNFaGwQwxnwJeAlnqcEHrLXveVO6iIiIyNkxPsP4uTmMnZlNxWsHWffifv70rTVMuWQUs5cVkZga63WJMsg8W+d7MGipQREREQknLY0drP3LXt5bUY3Pb5hyaR4zFxcqhA8BYb3O92BR+BYREZFwVH+0hXXP72PHmiNOCL8kj5lXKYRHMoVvFL5FREQkvNUfbWHdC/vYsVohPNIpfKPwLSIiIpGhawj3RxmmLShg1uJCbVcfQRS+UfgWERGRyFJ/pIU1z+5hV/lRYhOimLm4kGkLCoiO8XtdmpyBwjcK3yIiIhKZjh08yepn9rB/cy0JKTGULiti0sWj8Ed5tUq0nInCNwrfIiIiEtmqK+tZ9fRuaiobSEqPZebiQibNH0WURsLDjsI3Ct8iIiIS+ay1HNx6gvIX9lFT2UB8cjTTFxYw9bJ8YuK92rJFTne24Vs/MREREZEwZoyhcHImhZMzqd5Vz7oX97Hq6T1sePkAUy/PZ/qCAuKSdGJmpFD4FhEREYkQo0rSGFUyg6P7G1n/4n7KX9jHptcOMnPxaKYvLCA6VtNRwp2mnYiIiIhEqNrqJlb/eQ97K4fRl5MAABAvSURBVI6TkBJD2TVjmDg/F79fJ2YONs35RuFbREREhoea3Q2sfKqSmsoGUkfEM++6sYydlY0xxuvShg2FbxS+RUREZPiw1rJ/cy0rn97NiepmRoxOpuyaMYyekqkQPgh0wqWIiIjIMGKMoWhaFoVTMtm5+jBrn9vLcz/bxIiiFOZcM4bCyRkK4WFA4VtERERkCPH5DBMuzKVkzkh2rDpM+fP7+MtPKxg5JoWya8ZQOEkh3EuadiIiIiIyhAUDofdD+MkTbeQUp3LJR0sYMTrF69KGFM35RuFbRERE5JRgIMT2lTWseXYvLSc7mHxJHvOuKyYuUWuE9wfN+RYRERGR9/mjfEy+JI8LSkey9tm9bHqjit3rjnLhDWOZeGEuxqepKINBi0CKiIiIDCOx8VFcfHMJH72vjPTcBF7//Xae+P46ju5v9Lq0YUHTTkRERESGKWstO9cc4Z0nKmk92cHYmdmULisiKz/Z69IijqadiIiIiEivjDGMn5tD0bQsNry8n82vV7F7/THGTM+idFmRTsocABr5FhEREREA2ls62fR6FRWvHqS9JUDh5EzKri4ipzjV69LCnlY7QeFbREREpC86WgNsfrOKja8cpK2pk/wJ6ZQuLWLUuDStEd4DhW8UvkVERETOR2d7kC1vHWLD8gO0NnaQe0EqpUuLKNBGPX9D4RuFbxEREZH+EOgIsvWdGja8vJ+munZGjE5m9tIixkzPUgh3KXyj8C0iIiLSn07tlrnuxX00Hm9jxOhk5l0/loIJGV6X5jmFbxS+RURERAZCKBhix+ojrPnLHppOtFMwMZ0Lr7+A7MLhu0ShwjcK3yIiIiIDKdAZZMubh1j3wn7amju5oHQEc68tJm1EgtelDTqt8y0iIiIiAyoq2s+MRYVMnD+KjcsPsPGVA+xZf4yJ83MpXVZEUnqc1yWGHYVvERERETkvsfFRzL22mCmX5VH+/D62rqhm+8rDTL5kFLOWjCYxNdbrEsOGpp2IiIiISL9qPN5K+Qv72L7yMH6/YeoV+cxcXEh8UozXpQ0YzflG4VtERETES/VHWlj73F52rj1CdIyf6YsKmHllITFxQ2/yhcI3Ct8iIiIi4aC2uom1z+5l94ZjxCdHU3b1GCZdMgq/3+d1af1G4RuFbxEREZFwcnhvAyuf3E31rnpSR8Qz77qxjJ2VPSQ26lH4RuFbREREJNxYa9m/uZZ3n9pNXU0zI8ekcNENYxlVku51aedFSw2KiIiISNgxxlA0LYvCyRlsX3WYNc/s4akfbCB/QjpzrhlD7gVpXpc4oDTyLSIiIiKe6ewI8t5bh1j/0n5aT3Y6IfxDxeSOTfW6tHOiaScofIuIiIhEis72IFveOsSGl50QXjApgznXjCGnODJCuMI3Ct8iIiIikaaz3dmyfsNyJ4QXTc1k7nVjycpP8rq0Xil8o/AtIiIiEqk624Nsev0gG14+QHtrgJLSkcy9dgyp2Qlel9YtnXApIiIiIhErOtbP7CVFTL4kjw0vH2DTawfZve4oky4eRenVRRG7Zb3Ct4iIiIiErbjEaC68fizTFuRT/vw+tr5dzfaVNUy9Ip9Zi0cTlxTtdYnnRNNORERERCRiNBxrYc2zzpb1MbF+ZlxZyPSFBZ5vWX+200482dPTGDPDGLPKGLPRGFNujJnjthtjzH8ZYyqNMZuMMbO6POY2Y8wu93KbF3WLiIiIiLdSsxO48o7J3PLPc8gbn86aZ/fy+39eycZXDhDoCHpd3hl5MvJtjHkZ+JG19gVjzDLgbmvt5e7tLwPLgLnAj621c40xGUA5UApYYB0w21pb19vraORbREREZGg7sreR1c/s5uC2OhJTY7jiExMZPSVz0OsI65FvnACd4t5OBard29cBv7OOVUCaMSYXuApYbq094Qbu5cCSwS5aRERERMLLyDEpXPuVmXz4azNJyYonISXG65J65dXkmK8CLxlj/gPnF4CL3PY84GCX46rctp7aRURERETIG5/ODf842+syzmjAwrcx5hUgp5u77gMWAl+z1j5hjLkZ+DWwCDDdHG97ae/ude8E7gQoLCzsQ+UiIiIiIgNjwMK3tXZRT/cZY34HfMX98jHgV+7tKqCgy6H5OFNSqoDLT2t/o4fX/QXwC3DmfJ975SIiIiIiA8OrOd/VwGXu7QXALvf2M8Ct7qon84AGa20N8BKw2BiTboxJBxa7bSIiIiIiEcOrOd+fAX5sjIkC2nCniQDP46x0Ugm0ALcDWGtPGGO+Dax1j/uWtfbE4JYsIiIiInJ+PAnf1toVwN/MiLfOuodf7OExDwAPDHBpIiIiIiIDxqtpJyIiIiIiw47Ct4iIiIjIIFH4FhEREREZJArfIiIiIiKDROFbRERERGSQKHyLiIiIiAwShW8RERERkUGi8C0iIiIiMkgUvkVEREREBonCt4iIiIjIIFH4FhEREREZJArfIiIiIiKDROFbRERERGSQGGut1zUMGGPMMWC/Ry+fBRz36LWHA/XvwFL/Diz178BS/w4s9e/AUv8OrIHs39HW2uwzHTSkw7eXjDHl1tpSr+sYqtS/A0v9O7DUvwNL/Tuw1L8DS/07sMKhfzXtRERERERkkCh8i4iIiIgMEoXvgfMLrwsY4tS/A0v9O7DUvwNL/Tuw1L8DS/07sDzvX835FhEREREZJBr5FhEREREZJArf/cwYs8QYs8MYU2mMucfreiKVMWafMWazMWajMabcbcswxiw3xuxyr9PddmOM+S+3zzcZY2Z5W334McY8YIw5aozZ0qXtnPvTGHObe/wuY8xtXnwv4aiH/v2mMeaQ+x7eaIxZ1uW+e93+3WGMuapLuz4/umGMKTDGvG6M2WaMec8Y8xW3Xe/hftBL/+o93A+MMXHGmDXGmAq3f/+v2z7GGLPafS8+YoyJcdtj3a8r3fuLujxXt/0+nPXSv781xuzt8v6d4bZ7//lgrdWlny6AH9gNFAMxQAUwyeu6IvEC7AOyTmv7d+Ae9/Y9wP3u7WXAC4AB5gGrva4/3C7ApcAsYEtf+xPIAPa41+nu7XSvv7dwuPTQv98Evt7NsZPcz4ZYYIz7meHX50ev/ZsLzHJvJwM73X7Ue3hg+1fv4f7pXwMkubejgdXu+/JR4Ba3/f8Bn3dvfwH4f+7tW4BHeut3r78/ry+99O9vgRu7Od7zzweNfPevOUCltXaPtbYDeBi4zuOahpLrgAfd2w8CH+7S/jvrWAWkGWNyvSgwXFlr3wJOnNZ8rv15FbDcWnvCWlsHLAeWDHz14a+H/u3JdcDD1tp2a+1eoBLns0OfHz2w1tZYa9e7t08C24A89B7uF730b0/0Hj4H7vuwyf0y2r1YYAHwuNt++vv31Pv6cWChMcbQc78Pa730b088/3xQ+O5fecDBLl9X0fsHmPTMAi8bY9YZY+5020Zaa2vA+ccCGOG2q9/75lz7U/187r7k/lnzgVNTIlD/nhf3T/AzcUa39B7uZ6f1L+g93C+MMX5jzEbgKE6o2w3UW2sD7iFd++r9fnTvbwAyUf/26PT+tdaeev9+x33//sgYE+u2ef7+VfjuX6abNi0n0zfzrbWzgKXAF40xl/ZyrPq9f/XUn+rnc/NzYCwwA6gBfuC2q3/7yBiTBDwBfNVa29jbod20qY/PoJv+1Xu4n1hrg9baGUA+zmj1xO4Oc6/Vv+fo9P41xkwB7gUmAGU4U0m+4R7uef8qfPevKqCgy9f5QLVHtUQ0a221e30UeArnw+rIqekk7vVR93D1e9+ca3+qn8+BtfaI+w9CCPglH/x5WP3bB8aYaJxg+Adr7ZNus97D/aS7/tV7uP9Za+uBN3DmGqcZY6Lcu7r21fv96N6fijOtTf17Bl36d4k7ncpaa9uB3xBG71+F7/61Fihxz2COwTlR4hmPa4o4xphEY0zyqdvAYmALTl+eOvv4NuDP7u1ngFvdM5jnAQ2n/hQtvTrX/nwJWGyMSXf//LzYbZNunHbewfU472Fw+vcWd0WDMUAJsAZ9fvTIne/6a2CbtfaHXe7Se7gf9NS/eg/3D2NMtjEmzb0dDyzCmVf/OnCje9jp799T7+sbgdesc0ZgT/0+rPXQv9u7/GJucObTd33/evv5MBBncQ7nC85ZtDtx5nPd53U9kXjBOVO+wr28d6ofcea8vQrscq8z3HYD/Mzt881AqdffQ7hdgD/h/Nm4E+e3+0/1pT+BO3BO8qkEbvf6+wqXSw/9+3u3/zbhfNjndjn+Prd/dwBLu7Tr86P7/r0Y58+/m4CN7mWZ3sMD3r96D/dP/04DNrj9uAX4F7e9GCc8VwKPAbFue5z7daV7f/GZ+n04X3rp39fc9+8W4CE+WBHF888H7XApIiIiIjJINO1ERERERGSQKHyLiIiIiAwShW8RERERkUGi8C0iIiIiMkgUvkVEREREBonCt4iIiIjIIFH4FhGRs2KMmWqMOexu3SwiIn0QdeZDRERkqDDGJOFs4tGbWtv9JhD/BFwE/Cvwd/1dm4jIcKBNdkREhhFjzC+BT5/hsGxr7fHBqEdEZLhR+BYRGUaMMbOAN4H9OFuBt3RzWE8j3yIicp4051tEJMIYY/KNMR/ty2OtteuBG4BxwM+AOmvt8dMu1hjzhjGm6LTX/ZYxZrMxZqcx5s7T7vub40VE5G8pfIuIRJ6FwKy+Pthauxy4Hbga+PnZPMYYcxUwE5gBfAT4cF9fX0RkOFP4FhGJIMaYi4EfAjcaYzYaY8b05XmstX8AvgF8xhjzL2fxkGuB3wLRwJeAJ/ryusaYaX2tWURkKNBqJyIiEcRau8IYsxb4urV2y6l2Y8zbQHI3D/m6tfaVHp7r+8aYXOD/GmOqrLUP9PLSs4G1QC2wD/haH7+FeOAhY8x11tq9fXwOEZGIpRMuRUQijDFmLzDOWtvZD8+VA+wG9gLTrLUht/0N4JPW2n3GGB9wwFqbb4yJB/4H2Gmt/dcuz/P+8V3a/h64p5uXzQV2W2vnnG/9IiKRRiPfIiIRxBiTCTScHrz7MvJtjEkAngGagWtPBe9ujAd2AVhrW40x7wA57nN8y1rb7bQVa+1DwEOnvWah+5p9HTkXEYloCt8iIpFlDFB9eqO19pJzeRJ3NPsPwFTgCmvtnl4OnwnEGmP8OP9u/B1wlzFmBBBzLq+LE+S/YK199xwfJyIyJOiESxGRyLIdyDLGbDHGXHQez/Mj4DrgE9baVWc4dgbOXO3dwDvAg9baCqAMWHcuL2qtXa7gLSLDmUa+RUQiiLW2CTivudLGmLuAu4C7rbWPn8VDZuKE9C2ntZcBD55PLSIiw43Ct4jIMGKMuQ5n1PsF4FljzIRuDjtw2tcTcEbcTzdGK5aIiJwbhW8RkeHl2zhTDpe6l+5cgbOmdz2Atbbg9AOMMV8EKrs0vX+8iIj0TEsNioiIiIgMEp1wKSIiIiIySBS+RUREREQGicK3iIiIiMggUfgWERERERkkCt8iIiIiIoNE4VtEREREZJAofIuIiIiIDBKFbxERERGRQaLwLSIiIiIySP4/W76nB+2qpY0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize = (12,8))\n",
    "for i in range(10):\n",
    "    plt.plot(ts, coefs[:,i], label=label[i])\n",
    "\n",
    "plt.ylabel('coefficients')\n",
    "plt.xlabel('$t=\\sum|\\hat{\\\\beta}_j| \\\\rightarrow $')\n",
    "plt.title('Lasso paths - Sklearn')\n",
    "plt.legend()\n",
    "plt.axis('tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. LARS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Alg**: <br/>\n",
    "$initialization\\ \\hat{\\mu} = 0$ <br/>\n",
    "$repeat\\ t = 1,2,...,m:$<br/>\n",
    "$\\quad Cal\\ current\\ correlations\\ vector\\ \\hat{c}$ <br/>\n",
    "$\\quad Cal\\ X_{\\cal A},A_{\\cal A},u_{\\cal A}\\ and\\ a$ <br/>\n",
    "$\\quad updates\\ \\hat{\\mu}_{\\cal A}: \\hat{\\mu}_{\\cal A_+} = \\hat{\\mu}_{\\cal A} + \\hat{\\gamma} u_{\\cal A},$ <br/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lars_reg(X, y):\n",
    "    m = X.shape[1]\n",
    "    mu_hat = np.zeros_like(y)\n",
    "    w = np.zeros_like(X[0])\n",
    "    \n",
    "    for i in range(m):\n",
    "        # get residual with lars estimate mu_hat\n",
    "        residual = y - mu_hat\n",
    "        \n",
    "        # cal current correlations vector\n",
    "        c_hat = (X.T).dot(residual)\n",
    "        C_hat = np.max(np.abs(c_hat))      # max correlations\n",
    "        \n",
    "        # active set\n",
    "        A = set(np.where(np.isclose(np.abs(c_hat), C_hat))[0])\n",
    " \n",
    "        X_A = X[:, sorted(list(A))]*np.sign(c_hat[sorted(list(A))])\n",
    "        G_A = (X_A.T).dot(X_A)\n",
    "        G_A_inv = np.linalg.inv(G_A)\n",
    "    \n",
    "        A_A = (np.ones((1, len(A))).dot(G_A_inv).dot(np.ones(len(A))))**(-.5)\n",
    "        \n",
    "        w_A = A_A*G_A_inv.dot(np.ones(len(A)))\n",
    "        u_A = X_A.dot(w_A)   # equiangular vector\n",
    "        \n",
    "        a = (X.T).dot(u_A)\n",
    "        \n",
    "        # start cal step length: gamma_hat\n",
    "        tmp = [(C_hat-c_hat[j])/(A_A-a[j])  for j in range(m) if j not in A]\n",
    "        tmp.extend([(C_hat+c_hat[j])/(A_A+a[j])  for j in range(m) if j not in A])\n",
    "        \n",
    "        gamma_hat = C_hat/A_A if len(tmp)==0 else np.min([e for e in tmp if e > 0])\n",
    "        \n",
    "        mu_hat += gamma_hat*u_A\n",
    "        w[sorted(list(A))] += gamma_hat*w_A*np.sign(c_hat[sorted(list(A))])\n",
    "        \n",
    "        if len(tmp) == 0: break\n",
    "    return w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ -10.0098663 , -239.81564367,  519.84592005,  324.3846455 ,\n",
       "       -792.17563855,  476.73902101,  101.04326794,  177.06323767,\n",
       "        751.27369956,   67.62669218])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lars_reg(X, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 参考：\n",
    "- Forward-stagewise regression: LeastAngle_2002.pdf 或者 'The Elements of Statistical Learning' p60\n",
    "- LARS: LeastAngle_2002.pdf 及 'The Elements of Statistical Learning' p73"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
